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ABSTRACT 


The  algorithms  necessary  for  the  computer  aided  desian 
of  digital  filters  from  lowpass  prototypes  have  been 
developed  and  compared,  then  im.plemented  in  a  computer 
program  to  aid  in  the  instruction  of  digital  filter  design. 
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INTRODUCTION 


nr 


A  . 


'The  procedure  for  the  design  of  lowpass,  highpass 


bandpass 

and  bandstop 

digital  filters 

f  ror^ 

lowpass 

prototypes 

■i  c  we’’  ''  pc+-aV>T-; 

sri'^  i. s 

the  pret 

response  to 

WV  O  0  ^ 

desired 

ch. ar 2 c t e r i c  t i cs  .  This  transf  orr^tion  can  be  accoyrolished  in 
either  the  analog  (s^  dorain  or  the  digital  (z'  doir^ain 
(Figure  1 '  . 


i lowpass 


a  n  a  1  o  o  1= 
t  prototype : 
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ANALOG  TRANSFORMATIONS' 
(s  dor.ain) 


I  analog  | 
filter  i 

I _ i 


JL 


L _ ! 


BTLINFAR 
'  TRANSFORMATION 

i _ 

ii 

_ y _ 


!  BILINEAR  t 
! TRANSFORMATION i 
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f  z  d  o a  i  n  ' 
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I  1 


1  own  ass 


1 


'  final 

i  digital  (DIGITAL  TRANSFORMATIONS )  =:^  digital 

I  prototype  • 


I  e  ^  1  ^ 


liter 


Figure  1.  Two  irethods  for  digital  filter  design. 
Note  that  each  of  the  horizontal  arrows  actually 
represent  four  different  tr ansf orrat ions : 

1'  Icvpass  prototype  to  lowpass  filter 
2)  lowr>aES  nrototvpe  to  hichpass  filter 
3'  lowpass  prototype  to  bandpass  filter 
4'  lowpass  prototype  to  bandstop  filter 


1 


Therefore.  the  two  classical  desian  techniaues ,  analoa- 
dioital  and  digital-digital  design,  are  based  on  doing  the 
frequency  response  transformations  in  their  respective  anaioo 
or  digital  domains.  Since  the  theory  for  the  standard  filter 
types  ( Butterworth ,  Chebyshev,  and  elliptic)  is  developed  in 
the  analog  or  's'  domain,  the  two  methods  for  designing 
digital  filters  differ  basically  in  the  sequence  of  the 
design  steps.  In  the  analog-digital  method  the  analog 
prototype  is  transformed  and  then  digitized,  while  in  the 


digital-digital  m.ethod  the  analog  prototype  is  digitized  and 


then  transformed 


L"  „  ,  t  !  ■ 


Each  of  the  transf oPm.ations  (arrows  in  Figure  1)  involves 
substituting  a  function  of  's'  or  'z'  for  the  appropriate 


variap. 


in  the  filter  transfer  function.  The  purpose  oi 


this  thesis  is  tc  develop  efficient  algorithms  to  carry  out 
these  transformations  and  then  to  imple.ment  these  algorithm.s 
in  a  computer  program,  to  aid  in  the  design  of  digital 


filters . 


Computer-aided  design  of  digital  filters  offers  a 
significant  advantage  since  it  allows  the  digital  filter 
coefficients  to  be  calculated  quickly  and  accurately.  The 
accuracy  of  the  filter  coefficients  is  an  especially 
significant  point  since  the  frequency  response  of  the  final 
filter  is  very  sensitive  to  small  inaccuracies  in  the 
calculated  coefficients.  As  an  exam.ple,  compare  the 


freauency  response  for  the  filters  in  Figure  2  and  3. 


Figure 


2  shows  the  frequency  response  magnitude  for  a  1/2-dB  ripple 
Chebyshev  bandpass  filter,  and  Figure  3  shows  the  frequency 
response  when  one  of  the  coefficients  has  been  changed  by 
1/lOOOth.  _ _ 


i  • 

' 

i  . 

: 

I . «... 

1 

i  •  * 

1  ... 

j 

i 

• 

•- 

X 

• 

[\/\7 

. 

1  . 

' 

. r  •**““ 

i 

j 

; 

1 

/ 

1 

1 

I 

i 

i 

j 

/ 

/ 

i 

I 

1  1 

i  i 

i  f 

•.00  O.OO  •.!•  0,10  0,00  OJtft  •.90  0.90  •.«•  •.«• 

DIGITAL  FREQUENCY  (P/20XHZ) 


Figure  2.  Chebyshev  Bandpass  Filter 
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Figure  3.  Same  Filter  With  Denominator  Coefficient  of 

z  changed  to  -.67594 


B.  GENERAL  DISCUSSION 


Each 
of  two 
function 
namely 


of  the  transformations  involves  substitutina  a  ratio 
polynomials  of  s  or  2  into  the  filter  transfer 
which  is  also  a  ratio  of  two  polynomials  in  s  or  z , 


Y{y) 

H(y)=  - 

U  (y) 


p  (X) 

y=  - 

g  (x) 


(1.1) 


The  variables  'x'  and  ‘y‘  can  be  either  's'  or  'z'  dependina 
on  the  transf orm.ation  to  be  performed.  Therefore,  we  have 
either 


H  (y)  = 


ao  +  at  y  +  a2  y2  +  ...  +  ao>  y“  I 

_ ^ _  i 

bn  +  bi  V  +  bav^  +  ...  +  br,  y"  I  p(x) 

iy=  - 


or 


a  (X) 


(1.2) 


p ( X )  ^p ( X '  )  2 

ao  +  at  -  +  32 -  + 


g  (X ) 


(g (x) ) ^ 


H  (X)  = 


p  ( X )  ( p  ( X  )  )  ‘ 

bo  +  bi  -  +  bt -  + 

a  (x)  (g  (x)  )  * 


+  aii 


(p (x)  ) ' 


( a  (X )  )  "> 


+bn 


(p (x) ) " 


(a (x)  ) ' 


{1.31 


In  general.  the  order  of  the  numerator  'm',  is  less  than  the 
order  of  the  denominator  'n'  (except  for  the  transformations 
in  the  digital  domain  where  the  result  of  the  bilinear 
transf  orm.ation  is  to  make  n=m)  .  We  can  generalize  the 


4 


results,  however,  by  assur.ina  that  the  order  or  the  nu.xerator 
polynomial  Y(x)  is  equal  to  the  order  of  the  denominator 
polynomial  U(x)  and  making  the  additional  coefficients 

am  4  1  .  am  «  2 . ar  equal  to  zero.  Then  m:Ultiplying  both  the 

numerator  and  denominator  of  the  above  equation  by  (g(x))'’ 
gives 

an(a(xi)"  -t-  aip(xWg(xM''->  +  .  .  .  +  an(p(xM“ 

H(x)=  -  (1.4) 

b(i  (  g  ( X )  ) '’  +  bi  p  ( X  )  ( G  ( X )  ) "  -  ‘  +  .  .  .  +  bp  ( p  ( X )  )  “ 

or.  in  a  different  form. 


i=n 

Z  ai  p(x)‘a(xi'''* 
i  =  0 

K  ( X '  =  -  (1.5' 

i=n 

r  bi  p  (>:)•'  a  (x;  ’  - 
i  =  C 


Now  with  the  r.u.meratcr  and  deno-t inator  of  the  same  form., 
we  can  use  the  sam.e  transf orm.ation  for  both  the  num.erator  and 
deno.m.inator  polyno.mials  of  the  filter  transfer  function, 
namely . 


a'  (x)=  an  ' ai  'x+  a?  'x^^. 


.  +  a  p  '  x"  ' 


(1.61 

i=n 

=  I  ai  p(x)‘g(x)'’'‘ 
i  =  0 


In  Equation  1.6.  the  a.  are  the  coefficients  of  the  old 
nu.merator  or  denominator  polynomials  and  ai  '  are  the  new. 
transformed  numerator  or  denom.inator  polynomial  coefficients. 


Since  the  new  coefficients  a,  '.  are  a  linear  combination  of 


the  old  coefficients  ai  . 


and  the  coefficients  of  p!x;  and 


glx'.  .  this  suaaests  the  use  of  a  matrix  transformation; 


where 


a '  =  aA  ;  n ) 


(1.7) 


a  is  the  coefficient  row  vector  for  the  old 
numerator  or  denominator  polyno.mial  of 
lenoth=n-t-l  . 


a  =  Can  a  •  32  .  .  .  a-  1 


a'  is  the  coefficient  row  vector  fox'  the  new 
numerator  or  denom.inator  polynomial  of 
lencth=n  '  -t-1  . 


a  '  =  r  a  '  (,  a  ’  1  a  ’ :  .  .  .  ar  ’  ] 


A  <' n )  is  the  transformation  matrix  with  dimensions 
( n  +  l )  X  ( n  ' +1 )  .  The  elem.ents  of  A(n)  depend  upon  the 
coefficients  of  p(x)  and  c(x).  Note  that  if  the 
rcws  and  cclum.ns  of  Air.)  are  num.bered  0  to  n'  .  then 

i=n 

a  =  i.  a'A(n)ii  (l.&i 

1  =  0 


(Where  Ain't  -  denotes  the  element  in  row  i. 
column  1  of  m.atrix  Am).) 

is  the  order  of  the  transfer  function  before  the 
transformation 


IS  the  order  of  the  transfer  function  after  the 
transform, ation .  Note  that  n’  equals  n  times  the 
hiahest  order  of  p(x)  and  a(x). 


For  exarrlt  when  n= 

1  and  n ' =2 .  p ( x ) 

and  glx) 

are  second 

ordc-r  and  Equation  1 

. 7  becomes , 

a  e  'Alim,. 

Alt'o'  Af2)r2l 

_  — i  i 

1 

=  1 aO '  al ' 

a2’  j 

.  A  :  '  '  ■ 

A  1  2  !  •  A  (  2  '  -  .  i 

L_ 

_i  . 

- - 

1 

(1  .  ?  1 

€ 


In  the  case  of  the  lowpass  diaital,  hiahpass  diaitai .  and 
bilinear  transformations,  the  two  functions  p(x)  and  g(x;  are 
first  order.  Professor  P.H.  Moose  of  the  Naval  Postgraduate 
School  of  Monterey  has  developed  the  matrix  A(n)  for  these 
two  cases  by  expanding  Equation  1.6  using  the  binomial 
theorem.  (Reference  1).  This  method  cannot  be  used  for 
transf orm.ations  involving  higher  order  polynom.ials  such  as 
the  bandpass  and  bandstop  transf orm.at ions  and  it  is  difficult 
to  develop  because  intricate  manipulation  of  the  binom.ial 
coefficients  are  involved.  The  binomial  expansion  approach, 
therefore,  was  not  pursued. 

A  general  method  has  been  found  that  develops  the 
required  transf orm.aticn  matrices  in  an  iterative  manner.  In 
other  words,  the  matrix  used  to  transform,  a  numerator  or 
denominator  filter  polynom.ial  of  order  n  is  developed  frcr 
the  .matrix  for  the  polynomial  of  order  n-1.  This  approar.n  is 
easy  to  implem.ent  in  a  computer  program,  and  can  be  used  for 
any  of  the  polynomial  substitutions  needed.  The  next  chapter 
IS  concerned  with  developing  the  transformation  m.atrices  for 
each  specific  transformation.  However,  as  an  introduction  to 
tr.e  miethod  used,  consider  the  generic  case  where  both  p(xi 
and  a'vX)  are  second-order  polyno.m.ials .  Since  all  of  the 
analog  and  dioital  transformations  involve  substitutions  with 
polynom.iais  of  second  order  or  less,  this  second-order 
aeneric  case  will  be  all  that  is  needed  to  develop  tne 
specific  transf  orm.ations  later. 


/ 


With  p(x)=ax2  +  bx  +  c  .  a(x)  =  dx^ +  ex  +  f  .  and 
substituting  into  Equation  1.6  for  the  first  order  case  of 
n=l  IrejneiTibering  that  with  p{x)  and  g(x)  of  second  order,  the 
resulting  polynomial  will  be  second  order) ,  we  get 

a'  (x)  =  ao  '  +  ai  'X  +  a2  ’x*  =  an  (dx^  +  ex  +  f )  +  ai  (ax*+bx  +  c) 

(1.10) 

Using  Equation  1.8  or  1.9  it  is  easy  to  pick  out  the  matrix 
coefficients  A(l)i i  (i=0.1:  j=C,l,2)  for  the  first  order 

generic  case.  The  coefficients  are, 

if  e  dl  (1,11) 

A  (  i  )  =  i  c  b  a  i 

i_  _J 

Now  with  r.=2  (n'=4)  substitutinc  into  Equation  1.6  gives, 

a'  (X)  =  an  ■  +  ai  'X  +  az  'X‘  +  as  ’  x-  +  a^  '  x"*  (1.12) 

=  ac  (dx‘+ex+f)^  +  a;  (ax^+bx+c)  (dx^+ex+f)  +a2  (ax*+bx+c)‘ 

and  the  matrix  A{2)  is; 


1  f2 

2ef 

2df+e2 

2de 

d2  j 

A  (  2  ,  =  i  f  c 

ec+bl 

cd+f a+eb 

db  +  ea 

da  i 

(1.13) 

i  c^ 

2bc 

2ac+b" 

2ab 

a  i 

L_ 

_J 

To  see  how  this  matrix  is  derived  from  the  previous  one . 
consider  the  first  row  of  the  matrices  in  Equations  1.11  and 
1.13.  For  both  matrices,  this  row  is  multiplied  by  the  ac 


coefficient  in  the  old  polynom.ial  . 


For  n  = 


we  have  ao  p  i  x 


and  we  pick  out  the  coefficients  of  the  powers  of  x  (x'  )  to 
pet  the  matrix  elements  A(n)o i .  For  n=2 ,  we  have  a.  (p(x) )* 
and  we  do  the  same  thina  only  now  we  have  the  additional 
factor  of  p(x).  The  matrix  elements  can  be  calculated  by 
hand  in  this  fashion  but  it  requires  a  lot  of  algebra.  A 
simpler  method  is  to  use  the  matrix  elements  of  the  previous 
matrix  and  pet  the  appropriate  elements  by  convolvino 
coefficients.  For  example,  consider  the  sequence  of  elements 
in  row  0  (the  top  row)  of  Ad),  namely 


* 


fed 


(1.14) 


To  pet  the  elemients  for  row  0  of  A(2)  .  we  need  to  multiply 
(convolve  coefficients)  by  another  p(x).  Therefore,  with 
rows  numbered  C  tc  n.  and  colur.ns  number  0  to  n',  we  have  the 


f ol lowinc . 


For  row  0  column  0, 

A(2)oc  =  d  e  f 

f 

For  row  0  column  1 . 


the  appropriate  matrix  element  is 

=  f2 

e  d  (1.15) 

the  appropriate  matrix  element  is 


A  (2) 


e 

f 


f 

e 


=  2fe 


(1.16) 


For  row  0  colu.T.n  2.  the  appropriate  matrix  element  is 


A  (  2  r.  :•  =  d 


c 

f 


e 

e 


=  2fd+e* 


(1. 
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For  row  0  column  3.  the  appropriate  matrix  element  is 


A(2)o3=  def  =  2de 

fed  (1.18) 

For  row  0  column  4.  the  appropriate  matrix  element  is 

A(2)n«=  def=d2 

fed  (1.19) 

This  convolution  of  the  coefficients  of  the  old  matrix  with 
the  coefficients  of  p{x)  can  be  carried  out  for  every  row  of 
the  old  matrix  to  obtain  all  the  rows  of  the  new  matrix 
except  the  last  row.  In  the  case  of  the  last  row,  we  are 
multiplyina  by  an  additional  a(x)  rather  than  an  additional 
p(x)  and  so  must  obtain  the  last  row  of  the  new  matrix  by 
convolvina  the  coefficients  of  the  last  row  of  the  old  matrix 
with  the  coefficients  of  g(x). 

In  sumiTiary  we  have . 

(i)  for  row  0  to  row  n-1: 

-row  i  of  A(n)  =  row  i  of  A(n-l)  *  coefficients  of  p(x) 

(1.20) 


(ii)  for  row  n: 

-row  r.  of  A>'n;  =  row  n-1  of  A(n-l)  *  coefficients  of  g(x) 

(1.21) 

where  *  indicates  linear  convolution  of  the  coefficients. 

To  see  how  this  method  of  convolving  coefficients  works, 
consider  multiplying  two  polynomials: 

fli>:)=f2!xif3;x)  (1.22' 

1  0 


Since  'x'  is  the  independent  variable,  no  matter  what  letter 


we  use,  the  polynomials  will  remain  the  same.  Substitutinc 
' z"  '  ’  for  'x'  in  the  above  equation  gives  us. 


fl  (z-  M  =  f2  (z-  Mf  3  (z-  M 


(1.23) 


and  we  now  have  the  familiar  form  where  we  can  either 
multiply  the  two  functions  of  'z'  or  convolve  the  sequences 
found  by  taking  the  inverse  z-transf orm.  For  f2(2"*  )  and 
f 3 (z~  M  we  find 


f 2 (z- 1  ) 

=  ao 

+  a.  2-  ‘ 

+  a2  2*  ^  +  . 

.  .  +  an  Z‘ 

f3 (z- >  ) 

=  bo 

1 

N 

+  bz  z“  *  . 

.  .  +  bm  Z" 

It  is  now  clear  that  we  can  find  the  coefficients  of  the 
product  fl (x)  by  convolving  the  sequence  ( an  i  with  the 
sequence  ( bm I  to  produce  the  sequence  ( co  I  where. 


1  an  i 

TO 

o 

TO 

II 

,  az 

,  aa 

i  bm  1 

=  1  bo  ,  bi  , 

.  bz 

.  bn 

and  with 

q  =  m  +  n 

i  Cq  1 

II 

n 

a 

.  Cz 

.  Cn 

,  an  i 

.  bm  I 

(1.25) 

.  Cn  1 


Then  by  takina  the  z-transforr  of  the  sequence  I Co  I  and 
substituting  'x'  for  ’ z* *  '  we  get. 

fl(x)  =  C(  +  C)  X  +  .  .  .  +  Cq  x'’  (1.26) 
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We  now  have  a  simple  straight-forward  method  of 
developing  the  necessary  transformation  matrices  for  digital 
filter  design  by  picking  out  the  first-order  transformation 
matrix  using  Equations  1.6  and  1.11,  and  then  using  the 
convolution  of  coefficients  method  iteratively  on  the  rows  of 
each  successive  transformation  matrix  with  Equations  1.20  and 
1.21  until  the  appropriate  matrix  A(n)  has  been  generated. 
The  next  chapter  will  develop  the  transformation  matrices  for 
each  of  the  following  transf  orm.ations : 
bilinear  transformation 

dioital  lowpass  prototype  to  digital  lowpass  filter 

dicital  lowpass  prototype  to  digital  highpass  filter 

digital  lowpass  prototype  to  digital  bandpass  filter 

dioital  lowpass  prototype  to  digital  bandstop  filter 

analcc  lowpass  prototype  to  analog  lowpass  filter 

analoc  lowpass  prototype  to  analog  highpass  filter 

analoc  lowpass  prototype  to  analog  bandpass  filter 

analog  lowpass  prototype  to  analog  bandstop  filter. 

We  will  see  that  some  of  these  transformations  are  sim.ilar 
and  consequently,  we  will  not  have  to  develop 
separate  transformation  m.atrrces  for  each  case. 


II.  DEVELOPMENT  OF  SPECIFIC  TRANSFORMATIONS 

A.  INTRODUCTION 

The  individual  transformations  are  discussed  in  the 
following  sections.  The  specific  form  of  the  transformations 
used  in  each  case  are  from  Reference  2.  Since  the  final 
result  will  be  to  develop  a  com.puter  program  to  assist  in  the 
calculation  of  digital  filter  coefficients,  confining  the 
specific  transformations  to  the  form  in  Reference  2,  will 
allow  the  com.puter  program,  to  be  used  in  conjunction  with 
Reference  2  as  an  aid  in  designing  digital  filters. 

B.  BILINEAR  TRANSFORMATION 

The  bilinear  transformation  may  be  used  to  transform  an 
analog  transfer  function.  H s ;  ,  to  a  digital  transfer 
function  K(2).  The  transf orm.ation  from  H(s)  to  H(2)  is 
described  by 

H (2^  =  H(2) I 

j  z-1  (2.1) 

j  s=  - 

2  +  1 

Refemnc  to  Equation  1.1,  p(x)=z-l  and  g(x)=2  +  l.  Since  both 
of  the  transf  orm.ation  polynomials  p(x)  and  g(x)  are  first 
order,  the  bilinear  transformation  will  not  change  the  order 
of  the  transfer  function.  Therefore,  for  an  analog  transfer 
function  of  order  'n',  the  transformation  matrix  X{n)  will  be 
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(n+li  by  (n+1).  From  Equation  1.6,  the  numerator  or 
denominator  polynomial  of  the  transfer  function  H(z)  is  given 
by 

i=n 

a’(z)=ao‘+a''z+...+  an*z"=  I  as(z-l)Mz  +  l)”"‘ 

i=0  (2.2) 

where  the  unprimed  coefficients,  ai ,  come  from  the  analog 
transfer  function's  numerator  or  denominator  polynom.ial.  In 
matrix  form,  with  n=l  and  referring  to  Equation  1.11,  the 
first  order  matrix  A(l)  for  the  bilinear  transformation  is 

!  1  1  i 

A(l)  =  I  (  (2.3) 

1-1  i! 

Usinc  the  convolution  of  coefficients  method  developed  in 
Chapter  1,  the  transf orm.ation  matrix  A{n+1)  can  be  developed 
from  the  matrix  A(n).  The  results  for  the  bilinear 
transf  orm.ation  matrices  Ad)  through  A{5)  are  given  in  Table 
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TABLE  1 

FIRST-ORDER  TO  FIFTH-ORDER  BILINEAR  TRANSFORMATION  MATRICES 


!-l  ll 

L_  _l 


A(2)  = 

1 

2 

1  i 

i 

-1 

0 

1 

1  i 

1 

1 

-2 

1 

1  i 

— 

_j 

A(3)  = 

1 

3 

3 

Ti 

i 

-1 

-1 

1 

1  i 

i 

1 

-1 

-1 

1 

1  i 

1 

-1 

3 

-3 

1 

1  i 

— 

_j 

A(4)  = 

1 

4 

6 

4 

1  i 

1 

} 

1-1 

1 

-2 

0 

2 

1 

1  1 

i 

1 

( 

1  1 

1 

r\ 

-2 

0 

1  1 

i 

i-1 

1 

2 

0 

-2 

1  1 

1 

I 

1  1 

-4 

6 

-4 

1 

li 

L_ 

_j 

A(5>  = 

— 

1 

5 

10 

10 

5 

n 

1 

-1 

-3 

-2 

2 

3 

1 

1 ! 

1 

i  1 
{ 

1 

-2 

-2 

1 

1 

1 

1 

1-1 

1 

1 

2 

-2 

-1 

1 

ll 

1 

-3 

2 

2 

-3 

1 

1  1 

1 

-1 

5 

-10 

10 

-5 

1 

ll 

l_  -J 


Referring  to  Table  1  and  keeping  in  mind  the  convolution 
of  coefficients  method  that  produced  the  matrices,  the 
following  relations  are  apparent  and  will  be  useful  in 
computer  implementation. 

A(n)ot  (the  first  row)  ~  ^  j  ^  j=0.1....,n  (2.4) 

where  (  “  )  indicates  the  binomial  coefficient 

A(n)to  (the  first  column)  =  (-1)*  for  i=0,l,...,n  (2.5) 

A(n)in  (the  last  column)  =  1  for  i=0,l . n  (2.6) 

A(n)(n-iii  =  (-H"'JA(n)ii  for  i=0 , 1 ...,  integer  (n/2 ) +1  (2.7) 

A(n)i(  =  A(n-l)ii  +  A(n-l))<i-i>  for  i=0 . 1 , . . . , n-1  (2.8) 

j=l , 2 , . . . , n-1 

Notice  that  the  values  of  '  the  matrix  elements  for  the 
bilinear  transformation  are  whole  numbers.  This  is  im.portant 
from  a  num.erical  accuracy  standpoint  since  whole  numbers  can 
be  expressed  precisely  in  floating  point  form..  This  means 
that  the  bilinear  transf  orm.ation  will  cause  a  loss  of 
accuracy  of  the  original  coefficients  only  because  of  the  n+1 
multiplications  and  additions  used  to  calculate  each  new 
diaital  coefficient. 

C.  LOWPASS  TO  LOWPASS  DIGITAL  TRANSFORMATION 

The  lowpass  to  lowpass  diaital  transformation  is  used  to 
transform,  a  lowpass  prototype  digital  filter  H(z)  into  a 
lowpass  diaital  filter  H(z).  The  polynomial  substitution 
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used  is 


H(z)=  Hp  (z)  I 

z-a 

z=  -  (2.9) 

1-az 

where  the  subscript  p  indicates  the  prototype  digital  filter 

transfer  function.  The  design  constant  a  is  determined 

[References  2  and  3]  by 

sin  (Gr  /2  -  Gc  ’  /2) 
a  =  - 

sinlGr /2  -  Gr " I2\  (2.10) 

with  Sr  =  the  prototype  critical  frequency  (usually  pi/2)  and 
Gr  '=  the  desired  critical  frequency  of  the  lowpass  filter. 
Again,  referring  to  Equation  1.1  with  p(x)=  z-a  and  g(z)=  1- 
az ,  the  lowpass  digital  transformation  matrix  will  not  change 
the  order  of  the  transfer  function  and  the  transformation 
matrix  A(n)  will  be  (n+1)  by  (n+D.  Using  Equations  1.6  and 

l. 11  as  before,  the  first  order  transformation  matrix  A(l) 
for  the  lowpass  digital  transformation  is: 

i  1  -a  i 

A(l)=i  i  (2.11) 

i-a  1  i 
L_  _i 

Once  again  using  the  convolution  of  coefficients  method,  the 
hiaher-order  transformation  matrixes  A(n)  are  developed  from 
the  lower-order  matrix  A(n-l).  The  results  are  shown  in 
Table  2  for  the  first-order  to  fifth-order  transformation 

m. atrices . 
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TABLE  2 

FIRST-ORDER  TO  FIFTH-ORDER  LOWPASS  DIGITAL 
TRANSFORMATION  MATRICES 


A(l)  = 

1 

-a  1 

1 

1 

i  -a 

1 

li 

L. 

_J 

A(2)  = 

1 

-2a 

a*  i 

1 

i-a 

j 

(l+a*  ) 

1 

-a  i 

1 

i  a* 

-2a 

1 

1  i 

i 

f— 

— . 

A(3)  = 

X 

-3a 

3a*  -a^ 

i 

1 

-a 

(l+2a* ) 

(-2a-a3  )  a* 

i 

1 

1 

ia2 

1 

(-2a-a-"  ) 

(l+2a* )  -a 

( 

i 

i 

i 

3a* 

-3a  1 

i 

1 

_j 

A  '  4:  '  - 

A 

-4a 

"9  ^ 

i 

1 

1  —  n 

'  1  *  3a* 

^  (-3a-3o3 ' 

( 3a*  *a^  ' 

-C3  1 

j 

i  2 

1 

(  _  ~  ^ 

'  (-2a-2a*  ) 

a*  ' 

1 

i  _  ^  i 

1 

f  ^  2  ^  Ql 

( l+3a* - 

-a  j 

1 

1  ^  A 

-  4  o  3  6  a  * 

-4a 

1  1 

t _ 

Zj 

A  ?  = 

-5a 

lOa* 

-lOa* 

5o’ 

-0,5  1 

1 

f  —n 

f 1  +  ?a*  ' 

4-4a.-6a3  ) 

(6a*  +a^  ) 

1 -4a* -a*  ) 

a4  1 

j 

i  a‘ 

1 

(-2a-3a-'' 

)  {l  +  6a*+3oM 

( -3a-6a* -a^  ) 

(3a*  +2a^  ) 

-a*  1 

i 

1 

!  -a 

1 

(  3a*  +2a‘' 

)  ( -3a-6a* -a’ ) 

( l  +  6a*  +3a^  ) 

(-2a-3a3  ) 

1 

a2  I 

1 

1 

i  a'- 
j 

{  -4a*  -a-' 

)  (6a*  +4a^  ) 

( -4a-6a* ) 

(l+4a*  ) 

-a  1 

i 

i  -a' 

Ba'’ 

-lOa* 

10a* 

-5a 

1  ! 
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Referring  to  Table  2.  the  following 


relations  are 


observed  and  will  be  used  in  the  computer  implementation. 

A(n)oi  (the  first  row)  =  ^  ^  {-a)->  for  j  =  0,l,  ,n  (2.12) 

A(n)io  (the  first  column)  =  (-a)‘  for  i-0.1....,n  (2.13) 

A(n)it  =  A (n) ( n -I  . I r - .  .  for  i,j=o,l . n  (2.14) 

A(n)ii  =A(n-l)(i-i'<i-i'  (-a)A(n-l)(i-iii  (<^.15) 

for  i=l , 2 .  . n-1 

j  =  l  ,2 . n 


Note  that  the  elements  of  the  lowpass  digital 
transformation  matrices  are  not  whole  numbers  as  in  the 
bilinear  transformation  case.  but  are  polynomials  of  the 
design  constant  o.  Numerical  error  is  introduced  in  this 
case  by  both  the  n+1  additions  or  multiplications,  and  the 
error  in  calculating  the  matrix  elements  that  multiply  the 
original  coefficients. 


D.  LOWPASS  TO  HIGHPASS  DIGITAL  TRANSFORMATION 

The  lowpass  to  highpass  digital  transformation  is  used  tc 
transform  a  lowpass  prototype  to  a  highpass  filter.  The 
polynom.ial  substitution  for  this  transformation  is: 


H  (z) 


Hp  (z) 


z 


z-a 


1-az 


(2.16) 


Since  this  transf orm.ation  is  the  same  as  that  for  the  lowpass 
to  lowpass  digital  transf  orm.ation  except  for  the  minus  sign. 


we  can  make  the  substitution  of  '-z'  for  'z'  in  the  oriainal 


prototype  transfer  function 

and 

then 

use 

the 

sar  e 

transf  orm.at  ion 

m.atrices 

as  the 

lowpass 

to 

lowpass 

case . 

For 

the  lowpass  to 

highpass 

transformation , 

the 

design 

constant  a 

is  (References  2  and  3) 


cos  (Gc  /2  -  Gr  •  /2  ) 

a  =  -  (2.17) 

cos  (Gr  /2  +  Gr  ■  /  2) 


with  Gc  -  the  prototype  critical  frequency  (usually  pi/2)  and 
0c'  -  the  desired  critical  frequency  of  the  highpass  filter. 


E.  LCWPASS  TO  BANDSTOP  DIGITAL  TRANSFORMATION 

The  polynomial  substitution  which  transforms  a  digital 
lowpass  prototype  to  a  bandstcp  digital  filter  is  (References 
I  a d  3  : 


H  (z)  =  H-  (zM 


Iz  = 


2a  1-k 

z^  -  -  z  +  - 

1+k  1+k 

2a  1-k 

1  -  -  z  -  z  - 

1+k  1+k 


w  :  t  r. 


a  = 


cos  ( ©.  '  /  2  +  ©:  ' 


cos (Gu  ' /2  -  Gi ■ /2) 


ana 


K  =  tan(0c  /2)tan(0u  ’/2-e,  ' /2) 


(2.18' 


(2.19) 


(2.20 


with  once  aoaj.r. ,  Gr  =  the  prototype  critical  frequency 
(usually  pi '2)  .  O  '  =  the  desired  upper  critical  diaital 


frequency , 


and  Gt  ' 


=  the  desired  low  critical  diaital 

frequency . 

By  letting  b=  2a/(k  +  l)  and  c=  ( 1-k '  /  ( 1  +  k )  .  Equation  2.16  car. 
be  simplified  to 


H(z)  =  Hr,  (z)  I 

z^  -  bz  +  c 

z  =  - 

1  -  bz  +  cz* 


(2.21) 


Again.  using  Equations  1.6  and  1.11  the  first-order 
transformation  matrix  is 

A ( 1 )  =  1 1  -b  cl 

I  1  (2.22' 

•  —  —  K  ' 

I _  i 

*'cte  th  =  '  since  the  substitu'ri''’"  pel ync’'i al s  are  second 
order,  this  trar.sf ormatic"  will  double  the  order  of  the 

ra-rix  A '  r. '  will  be  (n.-*-!'  by  (2n-*-l'  .  The  higher  order 

transformation  matricies  are  again  developed  recursively  and 
the  results  for  the  first-order  to  third- order 

transformations  are  oiven  in  Table  3. 
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TABLE  3 

FIRST-ORDER  TO  THIRD-ORDER  BANDSTOP 
DIGITAL  TRANSFORMATION  MATRICES 


A(l)  = 

1 

-b  c 

c 

-b  1 

_ 

_i 

A(2)  = 

1 

-2b 

(b2+2c)  -2bc 

c*  i 
i 

c 

( -b-bc 

(l+b2+c2)  (-b-bc) 

1 

c  i 

i 

c2 

-2bc 

(b2+2c)  -2c 

1 

1  1 

L_ 

A{3)  = 

1 

-3c 

(3b2  +3c) 

(-b*-6bc)  ...  i 

1 

c 

(-2bc-b) 

(l+2b2  +bc  +2c, 

1 

{-b^ -2b-2bc-2bc*  )  ...  1 

1 

i  c^ 

1 

(-2bc-bc2)  (c3 +2b2  c+2c+b2  ) 

i 

c^ 

-3c2  b 

(3b*  c+3c*  ) 

1 

(Redundant  entries  not  shown  for  A(3)  due  to  space 
iir.itations .  Entries  not  shown  can  be  found  using  Equation 
2.25) 

Notice  the  increased  amount  of  calculation  needed  to 
determine  the  matrix  elements  for  the  bandstop  case.  Each 
ir.atrix  elerr.ent  is  a  function  of  'b'  and  'c'  which  are 
themselves  functions  of  the  design  constants  a  and  k.  Table 
3  only  gives  the  transformation  matrices  for  first-  to  third- 
order.  However.  to  illustrate  the  amount  of  calculation 
reauired.  the  matrix  element  A(5)25  is 

A(5)2-  =  (-12bc^ -6bc2 -8b3 c-6bc-6b^ -3b-3bc^ -6b3 c* -b* )  (2.23) 

Since  each  coefficient  in  the  final  filter  is  found  by 
multirlyinc  n  +  1  m.atrix  elements  by  the  n+1  original 
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coefficients  and  then  addina  the  n+1  products,  the  accuracy 
of  the  design  constants  and  the  original  coefficients  can  be 
critical  in  the  calculation  of  the  final  filter  coefficients. 


The  symmetry  relations  for  the  bandstop  transformation 

are : 


A(n)to  (the  first  colum.n)  =  c‘  for  i  =  0,l . n  (2.24) 

A(n)ii  =  A (n ) ( n -  1  1 ( 2 n - i 1  for  i=0,l,...,n  (2.25) 

j=0.1 . 2n 


Afn'ii  =  A(n-l))i  -  bA (n-l)i(i-)>  +  cA(n-l))((-2i  (2.26) 

for  i=0 , 1 . n-1 

j=0 , 1 , . . . , 2n 

F.  LOWPASS  TO  BANDPASS  DIGITAL  TRANSFORMATION 

Again,  the  polynom.ial  substitution  which  transforms  the 
lowpass  prototype  to  a  bandpass  filter  is 


H(z)  =  Hp  (2) 


2a)c  k-1 

-  -  z  +  - 

k+1  k+1 


2ak  k-1 

1  -  -  z  +  -  z* 

k+1  k+1 


(2.27) 


where  in  this  case  k  =  tan (0/2 ) cot (0r  ' /2  -  ©i  '/2)  and  a  is 
found  from,  the  same  expression  as  the  bandstop  case.  With 
the  simplification  of  b=2ak/(k+l)  and  c= ( k-1 ) / ( k+1 ) .  Equation 
2.27  becomes 
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K(2) 


Hp  (z' 


2^  -  bz  +  c 

2  =  -  -  (2.28) 

1  -  bz  +  cz* 

Since  this  is  the  same  form  as  the  bandstop 
transformation,  we  can  use  the  matrices  previously  developed 
if  '-z'  is  first  substituted  for  *2*  as  was  done  in  the 
lowpass  and  highpass  situation. 

G.  LOWPASS  TO  LOWPASS  ANALOG  TRANSFORMATION 

The  lowpass  to  lowpass  analog  transformation  is  used  to 
transform,  an  analog  lowpass  prototype  to  a  lowpass  analog 
filter.  The  polynomial  substitution  is  a  simple  one,  nam.ely 


H(s)  =  Hr  (s) 


(2.29) 


s  =  s/w. 


where  Wt  is  the  prewarped  critical  frequency  (wr =tan (0r /2 ) 

and  Or  is  the  desired  digital  critical  frequency 

(Gr =2  (pi ) fr  /fs  ) .  This  transformation  matrix  is  found  to  be 


A  ; ;  =  ]  C 


0  0 

1  /Wr  0 

0  1/Wr 


0 

0 

0 


1  /Wr  “ 


or  Ain  =  diacf  1,  l/w,  .l/wr^ 


,l/wc">  throughout 
(2.  30 
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Note  that  in  this  case,  each  new  filter  coefficient 
depends  only  on  multiplying  the  old  coefficient  of  s'  by 
(wi  .  Therefore,  any  loss  of  accuracy  due  to  the 
transformation  is  due  to  this  one  m.ultiplication  and  the 
error  introduced  in  calculation  of  the  matrix  element  (wc)'‘  . 

H.  LOWPASS  TO  HIGHPASS  ANALOG  TRANSFORMATION 

The  polynomial  substitution  for  the  lowpass  to  highpass 
transformation  is 

H(s)  =  Hp  (s) 

S=Wr/S  (2.31) 

Again  using  Equations  1.6  and  1.11,  the  first-order 
transformation  m.atrix  is 

Ti 

i  (2.321 

0  j 

_ I 

The  method  of  convolution  of  coefficients  can  be  used  to  get 
the  hioher-order  transformation  matrices.  However,  notice 
that  Equation  2.32  is  the  same  matrix  as  the  lowpass  to 
lowpass  transformation  with  the  reciprocal  taken  of  each  non 
zero  elem.ent,  and  then  each  row  reversed.  This  can  be  shown 
to  be  the  case  in  general .  so  that  the  lowpass  to  highpass 
transform.ation  matrix  is  given  by 
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0 


0 


(2.33) 


io 

I 

1 0  0 


V. 


1  i 

I 

r>  I 


A(n)  =  !  . 


1 0  wr "  - 1  .  .  .  0  0  i 

i 

Wr  "  0  ...  0  0  i 

l_  _l 


Again.  each  new  filter  coefficient  depends  only  on  one 
ir.ultiplication  although  the  calculation  of  the  matrix  element 
itself  will  require  calculations  involving  raising  the  design 
constant,  wr  .  to  the  (n-1)**'  power  where  i  is  the  power  of  s 
in  the  final  filter  transfer  function. 


I.  LOWPASS  TO  BANDPASS  ANALOG  TRANSFORMATION 

The  polynomial  substitution  for  transforming  a  lowpass 
analog  prototype  to  a  bandpass  analog  filter  is 

H(s)  =  Hp  (s)  ! 

i  s2  +  wo*  (2.34) 

|s  =  - 

Bs 

where  wo  is  defined  as  the  geometric  mean  of  the  passband; 

Wo  =  (WuWi  ) '^  ,  and  B  is  the  width  of  the  passband,  B=Wu-wi  . 

(wn  and  wi  are  the  upper  and  lower  critical  frequencies 
respectively).  Again  using  Equations  1.6  and  1.11  and 
picking  out  the  appropriate  matrix  elements,  the  first-order 
transformation  m.atrix  is  given  by 
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A(1 ) 


i  0  B 

Wo  ^  0 
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il 
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The  convolution  of  coefficients  technique  is  again  used  to 
generate  the  higher  order  transformation  matrices  with  the 
results  for  the  first-order  to  third-order  matrices  listed  in 
Table  4.  As  with  the  bandpass  digital  transformation,  the 
analoc  bandpass  transformation  will  double  the  order  of  the 
original  transfer  function  since  the  substitution  polynomial 
is  of  second  order. 

TABLE  4 

FIRST-ORDER  TO  THIRD-ORDER  BANDPASS 
ANALOG  TRANSFORMATION  MATRICES 
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J,  LOWFASS  TO  BANDSTOP  ANALOG  TRANSFORMATION 

For  the  analog  bandstop  transformation.  the  polynomial 
substitution  is  the  reciprocal  of  the  bandpass  case: 


H(s)  =  Hp  (s) 


Bs 


S*  +  We  * 


(2.36) 


Usino  the  same  method  as  before,  the  first-order 

transformation  matrix  is  found  to  be: 

[we^  0  I]  (2.37) 

A(l)  =  i  I 

io  B  0  1 

L_  _l 

Note  that  this  is  the  same  matrix  as  the  bandpass  case  but 
with  the  rows  in  reverse  order.  Using  the  convolution  of 
coefficients  technique  the  higher-order  matrices  are 
developed  and  the  results  for  the  first-order  to  third-order 
transformation  matrices  are  shown  in  Table  5.  Notice  that  the 
reversal  of  the  order  of  the  rows  compared  with  the  bandpass 
case  is  true  rn  general.  Therefore,  once  again  the  matrices 
are  over  half  empty  and  nonzero  elements  are  simple  functions 


t 


of  B  and  Wo  . 


TABLE  5 

FIRST-ORDER  TO  THIRD-ORDER  BANDSTOP 
ANALOG  TRANSFORMATION  MATRICES 


Ad) 


i  Wr.  2 


io 

i_ 


0 

B 


1  i 

I 

c  i 

_j 


A(2)  = 


Wo  * 
0 


0 

Bwc  -I 


2Wo2 

0 


0 

B 


0 


C  B2 


0 


0  I 

_J 


Wo  ^  0  2wc  ^ 


0  3wo  2  ^  0 


1 


A(3) 


to 

t 

I  0 


Bw..  “  0  2Bwc.  0 

0  B*  Wo  2  .0  B2 


B  0 

0  0  i 


io  0  0 

L_ 


B^ 


0 


K.  GENERAL  OBSERVATIONS 

Now  that  the  various  transf orr.ation  matrices  for  digital 
filter  design  have  been  developed,  it  is  possible  to  compare 
the  two  methods  of  diaital  filter  design  mentioned  earlier 
and  depicted  in  Figure  1.  For  this  comparison,  consider  an 
analoc  lowpass  prototype  filter  transfer  function  with 
numerator  of  order  m,  and  denominator  of  order  n  (m  less  than 
or  equal  to  n) . 

First  consider  the  digital  to  digital  direct  design 
method.  Transforming  the  analog  prototype  to  a  digital 
r'-otomvr.e  usino  the  b^^linear  transformation  will  require  m+1 
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multiplications  and  additions  per  digital  prototype 
coefficient  for  the  numerator  polynomial,  and  n+1 
multiplications  and  additions  per  digital  prototype 
coefficient  for  the  denominator  polynomial.  Then,  each  of 
the  digital  to  digital  transformations  will  require  2 (n+1) 
m.ultiplications  and  additions  since  a  result  of  the  bilinear 
transformation  is  to  make  n  =  m.  Therefore,  going  from  an 
analog  prototype  to  a  final  digital  filter  will  require 
(3n+m+4)  multiplications  and  additions.  If  lowpass  digital 
prototypes  are  developed  prior  to  the  design  process,  this 
can  be  reduced  to  2 (n+1)  since  then  only  the  digital  to 
dicital  transformations  will  be  used  in  designing  the  final 
filter  . 

Now  consider  the  analoo  desian  method  for  the  same  case 
of  an  r-  -order  numerator  and  n'  -order  denominator  lowpass 
prctctype  transfer  function.  For  the  lowpass  and  highpass 
case  each  new  coefficient  will  require  only  1  multiplication 
and  addition.  For  the  bandpass  and  bandstop  transformations 
a  simple  expression  for  the  number  of  calculations  per 
coefficient  is  difficult  to  find,  however  referring  to  Tables 
4  and  5.  the  number  of  multiplications  and  additions  per 
coefficient  is  always  less  than  half  the  order  of  the 
polynomial  to  be  transformed.  on  the  average.  This  gives 
approximately  (n-f-m)  /2  for  the  bandpass  to  bandstop  analoo 
transformations.  Finally  the  analog  filter  must  be  converted 
to  a  dicital  filter  using  the  bilinear  transformation  which 
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will  require 


(n+iri+2)  calculations  per  coefficient  for  the 


hiqhpass  and  lowpass  case  and  2{n+ir.+  l)  for  the  bandpass  and 
bandstop  case.  The  total  nur.ber  of  calculations  in  using  the 
analoq  design  method  then  is  (n+ir  +  4)  for  the  high/low  pass 
and  (  2 . 5  (n+ir. ) +2 )  for  the  bandpass/stop.  These  results  are 
summarized  in  Table  6. 


TABLE  6 

NUMBER  OF  MULTIPLICATIONS  AND  ADDITIONS  REQUIRED 
PER  COEFFICIENT  TO  FIND  THE  DIGITAL  FILTER 
COEFFICIENTS  FROM  A  LOWPASS  ANALOG  PROTOTYPE 
(WITH  AN  MT«  ORDER  NUMERATOR  AND  N^ "  ORDER  DENOMINATOR) 


Analog  Design 

high/low  pass  bandpass /stop 


analog  transformations  |  2  .5{n+m.) 

bilinear  transformation  !  n+ir,+  2  I  2(n+Tr)+2 


Total  I  n+m+4  2.5(n+ir.  )+2  j 

i  _ _ ! 


Digital 

Design 

high/low  pass 

bandpass /stop 

bilinear  transformation 
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i  1 

(  n+m+2ro  1 

1  ) 
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1  2(nel) 

i 
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1  2(n+l)  i 
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Referring  to  Table  6.  the  advantaae  of  one  method  over 
the  other  depends  upon  what  type  of  filter  is  being  designed 
and  the  order  of  the  prototype  transfer  functions  numerator 
and  denominator.  Consider  designing  a  fifth-order 
low/highpass  Chebyshev  or  Butterworth  filter.  In  this  case, 
the  order  of  the  numerator  polynomial  is  zero  and  the  analog 
method  will  have  a  slight  advantage  of  9  multiplications  and 
additions  versus  12  for  the  digital  design  with  precalculated 
prototypes  (and  19  for  digital  design  without  precalculated 
prototypes).  However,  if  a  bandpass /stop  elliptic  filter  was 
being  designed  with  n=5,  then  m=4  and  the  direct  digital 
design  method  (with  precalculated  prototypes)  has  an 
advantage  of  12  multiplications  and  additions  compared  to 
approximately  24  for  the  analog  method. 

A  final  consideration  in  favor  of  analog  design  is  the 
com.plfexity  of  the  elements  of  the  transf orm.ation  matrices. 
Since  the  elements  of  the  digital  transf  orm.ation  matrices 
require  much  more  calculation  than  the  elements  of  the  analog 
transformation  matrices  or  the  bilinear  transformation 
(com.pare  Tables  1-5),  the  problem  of  roundoff  error  is  much 
more  serious  in  the  case  of  the  digital  to  digital  design 
method . 

L.  IMPLEMENTATION 

The  final  step  is  to  implement  the  algorithms  developed 
previously  in  a  computer  program,  to  assist  in  the  teaching  of 
d-.cital  filter  design.  Since  the  program  is  intended  to  b- 
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used  for  instructional  purposes,  the  numerical  considerations 
discussed  earlier  are  not  considered  as  important  as 
conceptual  simplicity.  The  direct  design  m.ethod  was  chosen 
for  implementation  as  the  design  constants  can  be  found 
easily  without  reference  to  *’ prewarping The  completed 
proararr. .  called  'DFCADD'  (Digital  Filter  Computer  Aided 
Direct  Design) ,  is  contained  in  Appendix  A  with  its  own 
docum^ent  ati  on .  The  proarar  follows  the  procedures  set  forth 
in  Chapter  10  of  Reference  2.  Prototype  coefficients  art 
stored  within  the  program  for  first-order  to  fifth-order 
Butterworth  and  Chebyshev  {l/2,l.&2dB  ripples)  filters. 
While  second-order  to  sixth-order  elliptic  lowpass  prototype 
coefficients  are  calculated  within  the  program.  for  any 
combination  of  1/2, 1.2,  and  3  dB  passband  ripple  and  -20,- 
30,-40.-50  -60.  and  -10  dB  stopband  attenuation,  with  the 
help  of  a  prograr  developed  by  Professor  D.E.  Kirk  of  the 
Naval  Po.storaduate  School  and  References  4  and  5.  As  a  final 
option,  the  user  m.ay  enter  his/her  own  prototype  coefficients 
for  up  to  a  tenth-order  analog  prototype.  After  calculating 
the  digital  filter  coefficients  an  option  exists  for  the 
r'-oara-  to  create  a  data  file,  titled  '’filter”,  with  the 
dioital  filter  racnitude  and  phase  as  a  function  of  the 
dicital  frecuency  fror.  zero  to  pi.  This  data  file  can  then 
be  used  with  any  plcttinc  routine  to  generate  a  plot  of  the 
frecueniv  response  of  the  filter. 
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APPENDIX 

DFCADD  FORTRAN  PROGRAM 


FILE.  DFCADD  FORTRAN  A 


C 

C 


c 

c 


2 

1 

c 

c 

c 


9 

10 


KXKXKKKXXXXXXXXKKKVKXXKXXXXIfKXKXKXXXliXKKXXXXKXXKIOfliXXXRKiOfXXifXXK 
X  THIS  PROGRAM  ASSISTS  IN  THE  CALCULATION  OF  DIGITAL  FILTER  X 
X  COEFFICIENTS  ACCORDING  TO  THE  PROCEDURE  SET  FORTH  IN  CHAPTER  X 
X  10  (SUMMARIZED  OH  PG  <95)  OF  -  FIRST  PRINCIPLES  OF  DISCRETE  X 
K  SYSTEMS  AND  DIGITAL  SIGNAL  PROCESSING  -  BY  R.D.  STRUM  AND  X 
X  D.E.  KIRK  (REFERENCECZ)  )  X 

KXXXKXXXXXXXXXXKXXXXXXXXXXXXXKXXXXXXXXXXXXXKXXXXXXXXXXXXXKXKXKXX 
INTEGER  ORDER, TYPE, PASS, 0RDER2, END, AFORDF 
REAL  KPR,LAM,OMEGA(201,V(20),AO(2O),BO(2O),Bl(2O> 

COMPLEX  Z,NEWZ,ZM,ZD,HOFZN,HOFZD,HOFZ 
DOUBLE  PRECISION  TEMP,SUMR 

DOUBLE  PRECISION  R, SO , AOI , A02, BOI , B02, Bll , B12, HO 

DOUBLE  PRECISION  AOS, A04, BOS, BOA, BIS, BIA 

DOUBLE  PRECISION  SFREO, AFC, DFC, ALPHA, PI , P1D4, DFCD2 

DOUBLE  PRECISION  ALC , AUC, DLC, DUC, RK, B,C, DLCD2, DUCD2 

DOUBLE  PRECISION  APROT, BILINR, DPROT, LPHPTR, BPBSTR, DFLTR 

XXXKKXXXXXXXXXXXXXXXXXXXKXXXXKKKXHKKXXKXXXXXXXXKHXKKXKKKKXXKXKKX 

FORMAT  OF  ARRAYS i 

-APROT  CONTAINS  THE  ANALOG  PROTOTYPE 
APROT(NORD,PMR) 

-BILINR  IS  THE  ANALOG-DIGITAL  BILINEAR  TRANSFORMATION  MATRIX 
BILlNR(ORDER,ROH,COL) 

-DPROT  IS  THE  DIGITAL  PROTOTYPE  FROM  APROT  AND  BILINR 
DPROT(NORD,PWR) 

-LPHPTR  IS  THE  LPSHP  DIGITAL-DIGITAL  TRANFORMATION  MATRIX 
LPHPTR(ORDER,ROW,COL ) 

-DFLTR  IS  THE  FINAL  DIGITAL  FILTER 
DFLTR(HORD,PWR) 

-BPBSTR  IS  THE  BP8BS  DIGITAL-DIGITAL  TRANSFORMATION  MATRIX 
BPBSTR(ORD£R,ROH,COL) 

XXXXXXXXXXXXXXXXXXXKXXXXXXXXXXKXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXX 

DIMENSION  APROT(2,0-10) 

DIMENSION  Bill  HR (10,0. 10,0. 10) 

DIMENSION  DPRQT(2.0a0) 

DIMENSION  LPHPTRdO.OaO.OaO) 

DIMENSION  BPBSTR(10.0a0,0.20) 

DIMENSION  DFLTR(2,0:20) 

XXXXXXXXXXXXXXXXXXXXXKXKXXXXXXXXXXXXXXXXXXXXXKXXKXXXKXXXXXXXXXXX 
OPEN(NNIT=«, FI LE=’ FILTER  DATA • , STATUS^ ’UNKNOMN • ) 

XXX«XXX«XXXXXXXXXXXX»X>r»««>iXXXXXXXXXXKKKXXXXXXXKXXXXXXXXXXXXXXX 

PId.l4159265I5A98D+00 

SPI=S.14159 

PID4=PI/4.0D+00 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

CLEAR  APROT  PRIOR  TO  LOADING 
DO  1  K=l,2 
DO  2  L=0,10 
APROT(K,L)=0.0D»O0 
CONTINUE 
CONTINUE 

Xi<»»»X»*»««X«»««»»»»»«»»««X«»«XX«XXXXXXXXXXXXX»»X»»XXX«XXXXXXXXX 

xxxxxxxxxxxxoxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

BEGIN  INTERACTION 
CALL  EXCMS  ( 'CLRSCRN' ) 

HRITE(6.x)  'THIS  IS  A  PROGRAM  TO  CALCULATE  DIGITAL  FILTERS  USING’ 
WRJTE(6,«)  '  THE  DIRECT  DESIGN  METHOD.' 

I;RITE(6,») 

HRITE(6,X)  'XX»»»»X»*»X»X*«*»»*XXXXXXXXXXX»XXXXXX»XXXXXXXXXX' 

WRIT£(6,x)  'xHARNING-THIS  PROGRAM  IS  NOT  USER  FRIENDLY!!!  X' 
WRITE(6,X)  'X  YOU  MUST  ENTER  VALUES  OF  THE  APPROPRIATE  TYPEx' 
WRITE(6,x)  'x  (REAL  OR  INTEGER).  AND  PROPER  RANGE  x> 

ITRITE(8 , X)  'XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXX ’ 
HRITE(6,x) 

HRITE(6,X)  '  XX  ENTER  1  TO  CONTINUE  XX  • 

READ(5,X)  look 
CALL  EXCMS  ('CLRSCRN') 

HRITE(6,X)  'WHAT  TYPE  OF  FILTER  DO  YOU  HANTT’ 

IIRITE(6,X)  '  (1-4  ALLOW  ONLY  OP  TO  5TH  ORDER  PROTOTYPES)' 
WRITE(6,x)  •  l-BUTTERHORTH  • 

WRITE(6,»)  '  2-CHEBYSHEV  WITH  1^2  DB  RIPPLE  • 

WRITE(6,X)  '  S-CHEBYSMEV  WITH  I  DB  RIPPLE  ' 

HRITE(6,X)  •  4-CHEBYSHEV  WITH  2  DB  RIPPLE' 


NEHbuOlO 
NEW00020 
NEH00050 
NEM00040 
NEW00050 
NEH00060 
NEW00070 
NEW00080 
NEH00090 
NENOOIOO 
NEWOOllO 
NEM00120 
NEW00130 
NEH00140 
NEH00150 
NEWODUO 
NEH00170 
NEH00180 
NEH00I90 
NEW00200 
NEW00210 
NEW00220 
NEW00230 
NEW00240 
NEW00250 
NEW002B0 
NEW00270 
NEW00280 
NEW00290 
NEW00300 
NEW00310 
NEW00320 
NEW00330 
NEW00340 
NEW00350 
NEW0D560 
NEW00370 
NEMO  038  0 
HE)I00390 
NEW00400 
NEW00410 
NEW00420 
NEW00430 
NEW00440 
NEW00450 
NEW0C460 
NEII00470 
NEH00480 
NEM00490 
NEW00500 
NEWOOSIO 
NEW00520 
HEIT00550 
NE1T00540 
NEH00550 
NEW00S6O 
NEIT0057  0 
NEMO  0  58  0 
HEII00590 
NEW00600 
NEW00610 
NEW00620 
NEW00630 
NEW00640 
NEW00650 
NEW00860 
NEW00670 
NEH00680 
NEW00690 
NEH00700 
NEH00710 
NEH00720 
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FILEi  DFCADD  FORTRAN  A 


C 

lA 

15 


C 


WRITE(<,»)  ’  5-ElLIPTIC  (AllOWS  OP  TO  A  6TH  ORDER  PROTOTYPE)' 

HRITE(«,X)  •  -PASSBAND  RIPPIE  0 . 5, 1 . 0, 2 . 0, OR  3.0  DB' 

WRITE(<,»)  ’  -STOPBAND  ATTEMTUATION  -20,-30,-40,-50,-60,  • 

WR1TE(6,*)  '  OR  -70  DB  • 

WRITE(6,k)  '  6-OTHER  -  ALIOHS  YOU  TO  USE  UP  TO  A  lOTH  ORDER' 

HRITE(6,K)  ’  PROTOTYPE,  HOWEVER  YOU  MUST  ENTER  THE  ANALOG' 

WRITE(6,*)  '  PROTOTYPE  COEFFICIENTS  ' 

WRITE(6,*) 

WRITE(6,k)  '  »«  INTEGER  1-6  ««' 

READ(5,»)  TYPE 

KXKXXKKXXXKXKKKKKKKXKXXKXKKKXXKKXKXXKXKKKXXXXXXXXXXXXXKXXJfXKKKXK 
CALL  EXCMS  ('CLRSCRN') 

WRITE(6,X) 

WRITE(6,X)  'WHAT  IS  THE  ORDER  OF  THE  PROTOTYPE  YOU  WISH  TO  USET ' 
WRITE<6,X) 

WRITE(6,X)  •  -  MUST  BE  AN  INTEGER  IN  THE  RANGE  1-5  UNLESS  YOU' 
WRITE(6,X)  '  CHOSE  TO  ENTER  YOUR  OWN  PROTOTYPE  (THEN  1-10)' 
WRITE(6,X)  •  OR  YOU  CHOSE  AN  ELLIPTIC  FILTER  (THEN  1-6)' 
WRITE(6,X)  •  -  OR  ENTER  "0"  IF  YOU  WANT  HELP  SELECTING  THE  ORDER 
WRITE(6,X)  '  BY  CHANGING  SPECIFICATION  FREQUENCIES  TO  THE' 
WRITE(6,X)  •  CORRESPONDING  DIGITAL  PROTOTYPE  FREQUENCIES  ' 
WRITE(6,X)  '  NOTEi  YOU  WILL  MEED  PROTOTYPE  FREQUENCY  ' 

HRITE(6,X)  ■  RESPONSE  PLOTS  LIKE  THOSE  GIVEN  IN  R£F(2)  FOR 

WRITE(6,X)  '  BUTTERWORTH  FILTERS  (FIG  10. SI,  PGS  691B692)  ' 

WRITE{6,X)  ■  AND  CHEBYSHEV  FILTERS  (FIG  10.14,  PGS  696-698) 
READ(5,X)  ORDER 

IF((TYPE.GT.4) .AND. (TYPE.LT.8))THEN 
1F(0RDER.E0.1)THEN 

WRITE<6,X)  'ELLIPTIC  FILTERS  GIVEN  FOR  2ND-9TH  ORDER  ONLY' 
WRITE(6,X) 

GO  TO  10 
EHDIF 
ENDIF 


MUST  BE  AN  INTEGER  IN  THE  RANGE  1-5  UNLESS  YOU' 
CHOSE  TO  ENTER  YOUR  OUN  PROTOTYPE  (THEN  1-10)' 

OR  YOU  CHOSE  AN  ELLIPTIC  FILTER  (THEN  1-6)' 

OR  ENTER  "0"  IF  YOU  WANT  HELP  SELECTING  THE  ORDER' 
BY  CHANGING  SPECIFICATION  FREQUENCIES  TO  THE' 
CORRESPONDING  DIGITAL  PROTOTYPE  FREQUENCIES  ' 
NOTEi  YOU  WILL  MEED  PROTOTYPE  FREQUENCY  ' 
RESPONSE  PLOTS  LIKE  THOSE  GIVEN  IN  R£F(2)  FOR  ' 
BUTTERWORTH  FILTERS  (FIG  10. SI,  PGS  6918692)  ' 

AND  CHEBYSHEV  FILTERS  (FIG  10.14,  PGS  696-698)  ' 


NEH007SO 

NEW00740 

NEH00750 

NEH00760 

NEH00770 

NEH00780 

NEW00790 

NEW00800 

NEH00810 

NEW00820 

NEW008S0 

NEW00840 

NEW00850 

NEH00860 

NEHOOSTO 

NEW00880 

NEH00890 

NEW00900 

NEH00910 

NEH00920 

NEH009S0 

NEW00940 

NEW00950 

NEW00960 

NEW00970 

NEW0098U 

NEW00990 

NEWOlOOO 

NEWOlOlO 

NEW01020 

NEHOIOSO 

NEW01040 

NEWOlOSO 

NEW01060 


CALL  EXCMS  ('CLRSCRN') 

WRITE<6,X)  'DO  YOU  WANT  TO  ENTER  SPECIFICATION  FREQUENCIES' 
IIRITEC6,X)  'AS  ANALOG  OR  DIGITAL  FREQUENCIES?' 

HRITE(6,X)  '  1-ANALOG  ' 

HRnE(6,X)  '  2-DIGITAL  ' 

WRITEC6,X)  '  XX  integer  1  OR  2  xx  » 

READ(S,x)  AFORDF 
CALL  EXCMS  ('CLRSCRN') 

IFCAFORDF.EQ.DTHEN 

WRITE(6,X)  'WHAT  IS  THE  SAMPLING  FREQUENCY  IN  KHZ?' 
WRITE(6,X)  '  XX  DECIMAL  xx' 

READ(5,x)  SSFREQ 
SFREQ=DBLE<SSFREQ) 

ENDIF 

CALL  EXCMS  ('CLRSCRN') 

I'lRITECO.X) 

IJRITE(6.X)  'WHAT  TYPE  OF  PASSBAND  DUE  YOU  WANT?' 

I•!R^E(6,X)  '  l-LOWPASS' 

WRITE(6,X)  '  2-HIGHP.AS',' 

WRnE(6,x)  '  S-BANDFASS' 

WRITECA,*)  '  4-BANDSTOP' 

WRITE(6,X) 

WRITE(6,X)  '  XX  INTEGER  1-4  XX  ' 

READ<5,X)  PASS 

FORMULAS  FOR  HIGHPASS  AND  LOWPASS 
IFKFASS.EQ.l)  .  OR  .  ( PASS  .  EQ  .  2  )  )THEN 
CALL  EXCMS  ('CLRSCRN') 

IFtAFORDF.EQ.UTHEN 

WRITE(6,X)  'WHAT  IS  THE  CRITICAL  FREQUENCY  IN  KHZ?  ■ 
WRITE(6,X)  '  -  NEG  SDB  PT  FOR  BUTTERWORTH  ' 

WRITE(6,X)  '  -  RIPPLE  EDGE  FOR  ELLIPTIC  AND  CHEBYSHEV 

WRITE(6,X)  '  XX  DECIMAL  XX  ' 

READ(5,X)  SAFC 
AFC=DBLE(SAFC) 

DFC=(AFC/SFREQ)XPIX(2.0D+00) 

EHDIF 

IF(AFDRDF.EQ.2)THEN 

WRIT£(6,X)  'WHAT  IS  THE  CRITICAL  DIGITAL  FREQUENCY?’ 
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HRITE(6.«)  •  -  NEC  3DB  FOR  BUTTERWORTH  • 

HRITE(6,K)  '  -  RIPPLE  EDGE  FOR  ELLIPTIC  AND  CHEBYSHEV 

WRITE<6.I()  '  *»  DECIMAL  «ll  • 

READLS.K)  SDFC 
DFC'DBLECSDFC) 

ENDIF 

DFCD2*DFC/2.0D+00 

CALCULATE  THE  APPROPRIATE  VALUES  OF  THE  DESIGN  CONSTANT  ALPHA 
FROM  TABLE  10.7  REFERENCE  (2)  AND  REFERENCE  (3) 

IFCPASS.EQ. IJTHEN 

ALPHA-DSIN((PID4)-(DFCD2))/DSIN((PID4)+(DFCD2)> 

ENDIF 

IFCPASS.EQ. 2)THEN 

ALPHA=DCOS(CPID4)-(DFCD2))/DCOS((PIDA)+(DFCD2)) 

ENDIF 

KXKKAXKKVXXKXAXKKXXXXXXXXFXXXXIEXKXKXXKXKKKKKKKKKXKKXKKXXKKKKKMKK 

OPTION  TO  SHOH  CALCULATED  VALUE  OF  ALPHA 
CALL  EXCMS  ('CLRSCRN') 

WRITE(6,X)  'DO  YOU  WISH  TO  SEE  THE  CALCULATED  VALUE  OF  ALPHA!* 
WRITE(6,X)  •  XX  1-YES/2-NO  xx* 

READCS.Xl  LOOK 
IFCLOOK.EO. IJTHEN 
WRITE(6,x) 

WRITE<6,903)  ALPHA 
IIRITE(6,X) 

WRITE(6,X)  'XX  ENTER  I  TO  CONTINUE  XX* 

READCS.XJ  LOOK 
ENDIF 
ENDIF 

KXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXKXX 

BANDPASS  AND  BANDSTOP  FORMULAS 
IF((PASS.EQ.3) .OR. (PASS.EQ.AJJTHEN 
CALL  EXCMS  ('CLRSCRN') 

IFCAFORDF.EO. IJTHEN 

HRITEC6,XJ  'WHAT  IS  THE  LOWER  CRITICAL  FREQUENCY  IN  KHZT* 
WRITE(6,XJ  *  -  NEG  3DB  FOR  BUTTERWORTH* 

WRITE(6,*J  '  -  RIPPLE  EDGE  FOR  ELLIPTIC  AND  CHEBYSHEV  * 

WRITE(6,XJ  '  XX  DECIMAL  XX* 

READ(5,xJ  SALC 
ALC=DBLE(SALC) 

CALL  EXCMS  ( 'CLRSCRN') 

WRITE<6,KJ  'WHAT  IS  THE  UPPER  CRITICAL  FREQUENCY  IN  KHZ?* 
WRITE(6,X)  '  -  NEG  3DB  FOR  BUTTERWORTH* 

WRITE(«,XJ  '  -  RIPPLE  EDGE  FOR  ELLIPTIC  AND  CHEBYSHEV  * 

WRITE(6,X)  '  XX  DECIMAL  KK* 

READ(5,XJ  SAUC 
AUC=DBLE(SAUC) 

DLC=(ALC/SFREq)x(2.0D+00)»PI 
DUC=(AUC/SFREQ)x(2.0D+00 JXPI 
ENDIF 

IF(AFORDF.EQ.2)THEN 

WRITE(6,X)  'WHAT  IS  THE  LOWER  CRITICAL  DIGITAL  FREQUENCY!' 
WRITE(6,XJ  '  -  NEG  3DB  FOR  BUTTERWORTH  * 

URITE(6,X)  '  -  RIPPLE  EDGE  FOR  CHEBYSHEV  AND  ELLIPTIC 

WRITE(6,X)  '  XX  DECIMAL  xx  * 

READ(5,X)  SDLC 
DLC=DBLE(SDLCJ 
CALL  EXCMS  ( 'CLRSCRN* J 

WRITE(6,XJ  '  WHAT  IS  THE  UPPER  CRITICAL  DIGITAL  FREQUENCY!' 
WRITE(6,XJ  '  -  NEG  3DB  FOR  BUTTERWORTH  * 

WRITECS.X)  '  -  RIPPLE  EDGE  FOR  CHEBYSHEV  AND  ELLIPTIC* 

WRITE(6,X)  '  XX  decimal  XX  * 

READ(5,X)  SDUC 
DUC=DBLE(SDUCJ 
ENDIF 

DLCD2=DLC/2.0D+00 
DUCDZ=DUC/2.0D+00 
C  CALCULATE  APPROPRIATE  VALUES  OF  DESIGN  CONSTANTS  ALPHA  AND  K 
C  FROM  TABLE  10.7  REFERENCE  (2)  AND  REFERENCE  (3) 

ALPHA=DC0S(<DUCD2J+(DLCD2J J/DC0S((DUCD2J-(0LCD2J J 
IFIPASS.EQ.AJTHEN 

RK=DTAN(PIDAJXDTAN((DUCD2J-(DLCD2)J 
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C 


16 


B=((2.0D+00)kALPHA)/(RK+1.0D>00> 

C*((1.0D+00)-RK)/((1.0D+00)+RK) 

CNOIF 

IF(PASS.EQ.3)THEN 

RK=DTAN(PID4)/0TAN((DUC02)-(DLCD2)) 
B»((2.0D'»00)iiALPHAKRK)/(RK+(l  .0D«00)) 

C°(RK-(1 .0D+001)/(RK+(1  .OD-fOO)) 

ENDIF 

OPTION  TO  SHOH  CALCULATED  VALUE  OF  ALPHA  AND  K 
CALL  EXCHS  CCLRSCRN') 

HRITE(6,K)  'DO  YOU  HISH  TO  SEE  THE  CALCULATED  VALUE  OF  ALPHA  • 
HRITE(6,K)  •  AND  Kt  • 

NRITE(6,«)  >  KK  1-YES/2-N0  M* 

READ(5,X)  LOOK 
IFCLOOK.EO.llTHEN 
HRITE(6.«) 

HRITE(6,90S)  ALPHA 
HRITE(6,90<i)  RK 
HRITE(6,30 

HRITE(6.X)  '*»  ENTER  1  TO  CONTINUE  KK* 

READ(5,X)  LOOK 
ENDIF 
ENDIF 

XXXIIXKKXXXXKKXKItKKXXKKKKKXXXXKXXXKKXKKXXXXNKXKXXXXKXKXXXXKKXXXXII 

OPTION  TO  ASSIST  IN  DETERMINING  THE  PROTOTYPE  ORDER  BY  CHANOING 
SPECIFICATION  FREQUENCIES  TO  APPROPRIATE  PROTOTYPE  FREQUENCIES 
FROM  EQN  10.274  REFERENCEd) 

IFIORDER.EQ.OITNEN 
CALL  EXCMS  CCLRSCRN') 

IFCAFORDF.EQ.DTHEN 

WRITE<6,X)  'WHAT  IS  THE  ANALOG  SPECIFICATION  FREQUENCY  (KHZ)' 
HRITE<6,X)  'THAT  YOU  WANT  TO  CONVERT  TO  A  DIGITAL' 

WRITE(6,X)  'LOWPASS  PROTOTYPE  FREQUENCY?' 

WRITE(6,X)  '  XX  decimal  XX' 

READ<5,X)  F 

DF=(2.0)XSPIX(F/SSFREQ) 

ENDIF 

IF(AF0RDF.EQ.2)THEN 

WRITE(6,X)  'WHAT  IS  THE  DIGITAL  SPECIFICATION  FREQUENCY' 
WRITE<6,X)  'THAT  YOU  WANT  TO  CONVERT  TO  A  DIGITAL' 

HRITE<6,X)  'LOWPASS  PROTOTYPE  FREQUENCY?' 

WRITE(6,X)  '  XX  DECIMAL  XX' 

READ(5.X)  DF 
ENDIF 
X=COS(DF) 

Y=SINCDF) 

Z=CMPLX<X,Y) 

IFCPASS.EQ.DTHEN 

AS=REAL(ALPHA) 

HEWZ=(Z-AS)/((1 ,0)-(ASXZ)) 

ENDIF 

IF<PASS.EO,Z)THEN 

AS=REAL<ALPHA) 

NEWZ=(-1 .0)X((Z-AS)/((1 .O)-(ASXZ))) 

ENDIF 

IF(PASS.EQ.3)THEN 

BS=REAL(B) 

CS=REAL(C) 

NEWZ=(-1 .0)X(((ZXXZ)-(BSXZ)+CS)/((1.0)-(BSXZ)+(CSX(ZXXZ)))) 
ENDIF 

IF<PASS.E0.4)THEN 

BS=REAL(B) 

CS=REAL(C) 

NEWZ=((ZXXZ)-{BSXZ)+CS)/<(1 .0)-(BSXZ)+(CSX(ZXX2))) 

ENDIF 

PTHETA=ABS(ATAN(AIMAG(NEWZ)/REAL(NEHZ))) 

IFCPASS.EQ.DTHEN 

IFCF.GT.SAFOTHEN 

PTHETA=SPI-PTHETA 

ENDIF 

ENDIF 

IF(PASS.EQ.2)THEN 
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C 

C 

C 


C 


1F(F.LT.SAFC)THEN 

PTHETA'SPI-PTHETA 

ENDIF 

ENDIF 

IF(PASS.EQ.S)THEN 

IF( ( F . LT . SALC) . OR . ( F .GT . SAUC) )THEN 
PTHETA=SPI-PTHETA 
ENDIF 
ENDIF 

IF(PASS.EQ.4)THEN 

IFtCF.GT.SALO.AND.CF.LT.SAOCDTHEN 

PTHETA=SPI-PTHETA 

ENDIF 

ENDIF 

HRITE(6.905)  PTHETA 
MRITE(« ,K) 

HRITE(6,X)  'MANT  TO  CONVERT  ANOTHER  FREQUENCY!' 

MRITE(6,»)  •  »»  1-YES/2-N0  »«• 

READ(5,X)  LOOK 
IFCLOOK.EQ.DGO  TO  1« 

1F(L00K.EQ.2}THEN 

HRITEC6,X) 

WRITEC6,)<)  'WHAT  ORDER  PROTOTYPE  DO  YOU  WANT  TO  USE?' 
WRITE(6,)<)  '  »*  INTEGER  1-5  »»  ' 

READCS.x)  ORDER 
ENDIF 
ENDIF 

XKXXXXXKXXXKXXXXXKXXXXXXXXXXXXKKXXXXXXXXKIIKXXKKXXKXXKKKKKKKXKKKK 

LOAD  APRQT  MATRIX  WITH  APPROPRIATE  COEFFICIENTS 
APR0T(TYPE=1)  ARE  BUTTERWORTH  FILTERS 
IFITYPE.EO.DTHEN 
IF(0RDER.EQ.1)THEN 
APROT(1,0)=1.0D+00 
APROT(2,0)=l.flD+00 
APROT(2,U  =  1.0D+00 
ELSE1F<0RDER.E0.2)THEN 
APROT<l,0)=I.0D+00 
APROT(2,0)=1.0D+00 
APROT(2,1)  =  DS(3RT(2.0D+00) 

APROT(2,2)=l.OD+O0 
ElSEIF(ORDER.EQ.S)THEN 
APR0T(l,0)=l .OD+00 
APR0T(2,0)=1 .OD+00 
APR0T<2, 11=2.00*00 
APROT<2,2)=2.0D+00 
APROT(2,3)=l .OD+00 
ELSEIFCORDER.EQ.AjTHEN 
APROT(1,0)»1.0D+00 
APROT(2,0)'l,0D+00 
APROT(2,1)=2.61J126D+00 
APR0T(2,2)  =  3.Al<i21ADt00 
APROH2,3)=2.«13126D+00 
APROT(2,<i)  =  1.0D+0a 
ELSEIF(0RDER.EQ.5)THEN 
APROT(1,0)=1 .OD+00 
APROT<2,0)=1 .00+00 
APR0TI2,l)=3.236068D+00 
APROT( 2, 2) =5,2360680+00 
APR0T(2,3J=5.236068D+00 
APROT(2,‘i)  =  3.236068D+00 
APR0T(2,5)=I .OD+00 
ENDIF 
ENDIF 

APR0T<TVPE=2)  ARE  CHEBYSHEV  FILTERS  WITH  .5DB  RIPPLE 
IF(TYPE.EQ.2)THEN 
IFLORDER.EQ.DTHEN 
APROT(1,0)=2.862775D+00 
APROT(2,0)=2.862775D+00 
APR0T(2,1)=1 .OD+00 
ELSEIF(0RDER.EQ.2)THEN 
APROT(l,0)  =  l.<iS1388D+00 
APR0T(2,0)=1.516203D+00 
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APROT(2,l)  =  l.<i25625D+00 
APR0T{2,2)*X .00+00 
ELSEIF(0R0ER.Eq.3)THEN 
APROT(1.0)=0. 7156940+00 
APROT(2.0):0. 7156940+00 
APR0T<2,1)=1. 5348950+00 
APROTC 2, 2) -1.2529130+00 
APR0T(2, 31=1. 00+00 
ELSEIF(0RDER.EQ.4)THEN 
APROTd,  01=0.3578470+00 
APR0T(2, 01=0. 3790510+00 
APR0T(2, 1 1=1 . 0254550+00 
APR0T(2. 21=1. 7168660+00 
APR0T(2, 31=1. 1973860+00 
APR0T(2, 41=1 .00+00 
ELSEIF(OROER.E0.5)THEN 
APROTCl, 01=0. 1789230+00 
APR0T(2, 01=0. 1789230+00 
APR0T(2, 11=0. 7525180+00 
APROK  2, 21  =  1. 3095750+00 
APR0T(2, 31=1 .9373680+00 
APR0T(2, 41=1 .1724910+00 
APR0T(2, 51=1 .00+00 
EMDIF 
EHDIF 

C  APR0T(TYPE=3)  ARE  CHEBYSHEV  FILTERS  WITH  lOB  RIPPLE 
IF(TYP£.E«.3)THEN 
IFLORDER.EQ.llTHEN 
APROT< 1,01=1. 1965230+00 
APROTC 2, 01 =1.1965230+00 
APR0T(2, 11=1. 00+00 
ELSEIF(0R0ER.EQ.21THEN 
APROTC 1,01=0. 9826130+00 
APROTCa, 01=1. 1025100+00 
APROTC 2, 11=1. 0977340+00 
APR0TC2, 21=1. 00+00 
ELSEIFCOROER.E0.31THEN 
APROTCl, 01=0. 4913070+00 
APROTC2, 01=0. 4913070+00 
APR0TC2, 11=1. 2384090+00 
APROTC 2, 2 1=0. 9883410+00 
APR0TC2, 31=1 .00+00 
ELSEIFC0RDER.EQ.41THEN 
APROTCl, 01=0. 2456530+00 
APROTC 2, 01 =0.2756280+00 
APR0TC2, 11=0. 7426190+00 
APR0TC2, 21=1. 4539250+00 
APR0TC2, 31=0. 9528110+00 
APR0TC2, 41=1 .00+00 
ELSEIFCOROER.E0.51THEN 
APROTCl, 01=0. 1228270+00 
APROTC 2, 01 =0.1 22827 0+00 
APROTC2. 11=0. 5805340+00 
APROTC 2. 2 1=0. 974396 0+00 
APR0TC2, 31=1 .6888160+00 
APROTC 2, 4 1=0. 9368200+00 
APR0TC2, 51=1 .00+00 
EMDIF 
EMDIF 

C  APROTCTYPE=4)  ARE  CHEBYSHEV  FILTERS  HITH  2DB  RIPPLE 
IFCTYPE.EO.AITHEM 
IFCORDER.EO.llTHEH 
APROTCl, 01=1. 3075600+00 
APROTC 2, 01=1. 3075600+00 
APR0TC2, 11=1. 00+00 
ELSEIFC0RDER.E0.21THEN 
APROTC 1,0 1=0. 5058030+00 
APROTC 2, 01=0. 6367680+00 
APR0TC2, 11=0. 8038160+00 
APR0TC2, 21=1 .00+00 
ELSEIFC0RDER.EQ.31THEN 
APROTC 1,01=0. 3268900+00 
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APR0T(Z,0) =0.5268900+00 
APROTC 2, 1)=1. 0221900+00 
APR0T(2, 21=0.7578220+00 
APR0T(2, 51=1 .00+00 
ELSEIFC  OROER . EQ . AITHEN 
APR0T( 1,0 1=0.1659450+00 
APR0T(2, 01=0. 2057650+00 
APR0T(2, 11=0.5167980+00 
APR0T(2, 21=1. 256 4820+00 
APR0T(2, 51 =0.7162150+00 
APR0T(2, 41=1 .00+00 
ELSEIF(0R0ER.EQ.51THEN 
APROTC 1,01=0. 0817250+00 
APR0T(2, 01=0. 0817250+00 
APR0T( 2, 11=0. 4595490+00 
APR0T(2, 21=0. 6954770+00 
APR0T(2, 51=1. 4995450+00 
APROT( 2, 41=0. 7064610+00 
APR0T(2, 51=1. 00+00 
ENOIF 
ENDIF 

APR0T(TYPE=51  ARE  ELLIPTIC  FILTERS 

ELLIPTIC  COEFFICIENTS  ARE  CALCULATEO  USING  A  PROGRAM  OEVELOPEO 
BY  PROFESSOR  0.  E.  KIRK,  NAVAL  POSTGRAOUATE  SCHOOL,  MONTEREY, 
CALIFORNIA  95945 
A02=0. 00+00 
002=0.00+00 
012=0.00+00 
S0=0. 00+00 
IF<TYPE.E0.51THEN 

CALL  EXCMS  CCLRSCRN'l 

WRITE(6,*)  'WHAT  PASSBANO  RIPPLE  00  YOU  WANT?  (IN  001' 

HRITE(6,*1  '»**  MUST  BE  EITHER  0.5,  1.0,  2.0,  OR  5.0  KKK* 

READ<5,*1  ASUBP 

HRITE(6,X1 

MRITE(6,X) 

WRITE<6,>()  'WHAT  STOPBANO  ATTENUATION  00  YOU  WANTt  (IN  OBI* 
HRITE(6,*)  'K**  MUST  BE  -20,-50,-40,-50,-60,  OR  -70  »KK' 
WRITE<6,*1  '**K  00  NOT  INCLUOE  A  DECIMAL  UK*' 

READCS,*)  NATTN 
NATTN=-1»NATTN 
ASUBA=REAL(HATTtn 
IFCASUBP.EO.O.SITHEN 
IF<NATTN.EQ.201THEN 

IF<0R0ER.EQ.2)SR=2. 76261 
lF<ORDER.Eq.I}SR=l. 42189 
IF<0RDER.E<}.4)SR  =  1 .15188 
IF(0RDER.Eq.5)SR=l .04465 
lF(0RDER.Eg.61SR=l . 01555 
ENOIF 

IFLNATTN.EO.SOITHEN 

IFC0R0ER.EQ.2)SR=4. 80880 
IF(0P.DER.Eq.31SR  =  l.  92522 
IF(0RDER.E0.41SR=1 .52446 
IF(0RDER.Eq.5)SR=l .12912 
IF(0R0ER.Eq.61SR=l. 05594 
ENOIF 

IF(HATTN.Eq.401THEN 

IF( OROER. Eq. 2 )SR=5. 848925 
IF(0R0ER.Eq.51SR=2. 71147 
lF(0RDER.Eq.41SR=l .62842 
IF(0RDER.EQ.51Sfi=l .27264 
IF(0R0ER.Eq.6]SR=l .12697 
ENOIF 

IF(NATTN.Eq.50)THEN 

IF(OROER.Eq.21SR=I5. 06069 
IF{ORD£R.E0.51SR=3. 90450 
IF(0R0ER.Eq.4)SR=2. 06924 
IF( OROER. Eq.51SR=l .48469 
IF(0R0ER.Eq.61SR=l. 24101 
ENOIF 

IF(NATTN.Eq.601THEN 
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IF( ORDER. EQ.2)SR:26.7«2Sa 
IF( ORDER. EQ.S>SR*5.«7937 
1F( ORDER. EQ.4)SR:2. 68325 
IF( ORDER. Eq.5)SRxl. 77664 
IF(ORDER.EQ.6)SR<1.4013S 
ENDIF 

IF(NATTN.EQ.70)THEM 

IF(ORDER.Eq.2)SR-47. 58053 
IF(0RDER.Eq.3)SR°8. 30124 
IF(0RDER.E«.4)SR*3. 52135 
IF( ORDER. EQ.5>SR<2. 16377 
IF(0ROER.Eq.6)SR<l. 61413 
ENDXF 
ENDIF 

IF(ASUBP.EQ.1.0)THEN 

1F(NATTN.EQ.20)THEN 

lF(0RDER.Ee.2)SR:2. 32474 
IF( ORDER. Eg.3)SR-l. 30797 
IF(0RDER.EQ.4)SR:1 .09029 
IF(0RDER.E4.5)SR<1. 02826 
IF(ORDER.EQ.6)SR^l. 00902 
ENDIF 

IF(NATTN.EQ.30)THEN 

IF( ORDER. EQ.2)SR-4. 00423 
IF(ORDER.EQ.3)SR-l. 73254 
IF( ORDER. EO. 4 )SR-1. 25040 
IF(0RDER.Eq.5)SR=l .09554 
IF( ORDER. Eq.6)SR=l. 03799 
ENDIF 

IF(NATTN.Eq.40)THEN  , 

IF< ORDER. Eq.2)SR'7: 04488 
IF(0RDER.Eq.3)SR-2. 41619 
IF( ORDER. EO. 4 )SR=1 .51549 
IF(0RDER.Eq.5)SR=l. 21868 
IF(  ORDER.  E4. 6  )SR.=  1. 09887 
ENDIF 

IFtNATTN.EO.SOlTHEN 

IF< ORDER. Eq.2)SR=12. 48471 
IF(ORDER.£0.3)SR=3. 46061 
IF<0RDER.Eq.4)SR=l. 90819 
IF<0RDER.Eq.5)SR*l. 40723 
IF( ORDER . Eq .6  JSRsl .19891 
ENDIF 

IF(NATTN.Eq.60)THEN 

IF( ORDER. Eq.2>SR:22. 17688 
IF(0RDER.Eq.3)SR=5. 02121 
IF(0RDER.Eq.4)SR*2. 46079 
IF<0RDER.Eq.5)SR»l. 67161 
IF(0RDER.Eq.6)SR=1.34354 
ENDIF 

lF(NATTN.Eq.70)THEN 

IF(0RDER.Eq.2)SR=39.42285 
IF( ORDER. £0. 3 )SR=7. 33052 
IF(0RDER.Eq.4)SR=3. 21902 
IF(ORDER.Eq.5)SR=2.0Z567 
IFL0RDER.Eq.6)SR=l. 53842 
ENDIF 
ENDIF 

IF(ASUBP.Eq.2.0)THEN 

lF<NATTN.Eq.20)THEN 

IF(0RDER.Eq.2JSR=l. 94332 
IF(0RDER.Eq.3)SR=l. 20808 
IF( ORDER. Eq. 4 )SR=1 .05569 
IF10RDeR.Eq.5JSR=l .01567 
IF( ORDER. Eq. 6 )SR*1. 00447 
ENDIF 

IF(NATTN.Eq.50)THEN 

I F( ORDER. Eq.2)SR=3. 29235 
IF(DRDER.Eq.3)SR=l. 55690 
IF(0RDER.Eq.4)SR=l. 18280 
IF(0RDER.Eq.5)SR-l. 06594 
IFCORDER. Eq.67SR-l . 02460 
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ENDIF 

t  IF(NATTN.EQ.AO)THEN 

IF(0RDER.e0.2)SR*5.7<107 
IF( ORDER . EO . 5)SR=2 . 1392S 
IF( ORDER.  E«.<i)SR=1.408A2 
IF(0RDER.EQ.5)SR°1. 16811 
1F( ORDER. EQ.6)SR*1. 07316 

•  ENDIF 

1F(NATTN.E0.50)THEN 

IF(ORDER.E0.2)SR=10. 19181 
1F(0RDER.EQ.3)SR=3. 06137 
IF(0RDER.E8.9)SR»1. 75285 
lF(0RDER.Eq.5)SR<l. 33243 
lF(0RDER.Ea.6)SR:l. 15868 
ENDIF 

lF{NATTN.Efl.60)THEM 

1F(0R0ER.EQ.2)SR»18. 09398 
1F(ORO£R.E8.3)SR=4. 39729 
1F(0RDER.E«.4)SR=2. 24440 
1F( ORDER. EQ.S)SR-1. 56860 
1F(0RDER.EQ.6)SR=1. 28693 
ENDIF 

IF<NATTN.E0.70)THEN 

IFC  ORDER.  E<}.2)SR  =  32. 15951 
IFt  ORDER . E8 . 3)SR=6 . 40894 
IF(0RDER.E8.4)SR=2. 92363 
IF(0RDER.E8. 5)58=1.88906 
IF(0RDER.Ea.6)SR=l. 46330 
ENDIF 
ENDIF 

IF(ASUBP.E8.3.0)THEN 

IF(HATTN.E8.20)THEN 

IF(0RDER.E8.2)SR=I .73915 
IF(0RDER.E8.3)SR=1 .15516 
IF(0RDER.E8.4)SR=1 .03853 
IF<0RDER.E8.5)SR=1. 00996 
IF(QRDER.E8.6)SR=1. 00260 
ENDIF 

IF(NATTN.E8.30)THEN 

)  IF(0RDER.E8. 2)58=2.90352 

IF<0RDER.E8.3)SR=1. 45814 
IF(QRBER.E8.4)SR=1. 14542 
IFCORDER. £8.5)58=1.05020 
IF<0RDER.E8.6)SR»I. 01783 
ENDIF 

V  IF(HATTN.E8.40)THEN 

IF(ORDER. £8.2)58=5.05584 
IF<0RDER.E8. 3)58=1.98022 
IF(0RDER.E8.4)SR=1 .34663 
IF( ORDER. E8.5)SR=l .13934 
1 F( ORDER. E8.6)SR=1 .05892 
ENDIF 

IFCNATTN. EQ.50)THEN 

I F( ORDER. £8.2)58=8.93009 
IF( ORDER. £8. 3)SR=2. 79862 
1F(0RDER.E8.4)SR=1 .66148 
IF( ORDER. £8 .5)58=1 .28850 
I F( ORDER. E8.6)SR=1. 13534 
ENDIF 

IFCNATTN. E8.60)THEN 

IFCORDER.EO. 2)58=15.84610 
IFCORDER.E8.3)SR=4.03471 
IFCORDER.EO. 4)58=2.11596 
1FC0RDER.E8.5)SR=1 .50711 
1FC0RDER.E8.6)SR=1 .25325 
ENDIF 

"  IFCNATTN. E8.70)THEN 

IFC ORDER. E8.2)SR=28. 15950 
IFCORDER.EO. 3)58=5.87251 
IFCORDER.EO. 4)58=2.74749 
IFCORDER.EO. 5)58=1.80680 
IFCORDER.EO. 6)58=1.41799 
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17 


19 


ENDIF 

EHDIF 

SPI=J.1A159 
RK=1 .0/SR 
N=DRDER 
EN=N 

KPR=SQRT(1.0-RKkk2) 

Q0=0.5*(1.0-S«RT(KPR))/'(1.0+S<}RT(KPR)) 

Q-gO+2 .  OAQOUXS+IS .  0xg0i»t9+150 .  OXQOKKl! 

D=(10.D«x(a.lKASUBA)-l .0)/(10.0XX(0.1«ASUBP)-1.0> 

LAM=(1 .0/(2.  OXEN)  )XALOG((  10.  OXK(0.  OSKASUBPl-fl.  0)/(  10.  OXK(  0.05 
C  XASUBP)-1.0)) 

ISTOP'O 

DENSUM^O.O 

ENSUM°SINH(LAH) 

M=1 

EM=M 

ENUM=((-1.0)xxH)x(axx(Mx<M*l)))xSINH((2.0xEM+1.0)xiAM) 

DEN-2. Ox((-1.0)XXH)X(gxx(Mxx2))xCOSH(2.0XEnxLAM) 

ENSUM-ENSUM+ENUn 

OENSUM^DENSUN-fDEN 

IF(M.GE.2)THEN 

RN^ABSCENUM/ENSUM) 

RD=ABS(DEN/<1 .0+DENSUM)) 
IF((RN.LE.1.0E-8).AND.(RD.LE.1.0E-8))THEN 
IST0P=1 
ENDIF 
ENDIF 
M=M+1 

IFdSTOP.EO.OITHEN 
GO  TO  17 
ENDIF 

SIGMAO=ABS( <  2 . OX(QXXO . 25)XENSUM)/( 1 . 0+DEHSUM) ) 

W=SqRT( ( 1 . 0+RKX(SIGMA0xx2) )x( 1 . 0+ (SIGMA0XK2/RK) ) ) 

IN=M0D(N,2) 

IF<IN.EQ.O)THEN 

IR=N/2 

ELSE 

IR=<N-l)/2 

ENDIF 

DO  18  1:1. IR 
IFdN.EO.OJTHEN 
EMU=I-O.S 


ENDIF 

ISTOP^O 

DENSUMsO.O 

EHSUM=SIH(PIXEMU/EN) 

M’l 

EM=M 

ENUM=((-1 .0)xxM)x(qxx(Mx(M+l)))xsiN((2.0XEM+l .0)xPIxEMU/EN) 

DEN=2.0X((-1.0)XXM)X(CXX(MXX2))XCOS(2.0XEMXPIXEMU/EN) 

ENSUH=ENSUM+ENUM 

DENSUM=DENSUM+DEN 

IF(M.GE.2)THEN 

RN=ABS(ENUM/ENSUf.;) 

RD  =  ABS(DEN/d  .0  +  DENSUM)) 

IFdRN.LE.l  .0E-8).AND.(RD.LE.l.OE-8))THEN 
ISTOP=I 
ENDIF 
ENDIF 
«=M+1 

IFdSTOP.Eq.O)THEN 
GO  TO  19 
ENDIF 

OMEGAd  )  =  (2 . 0X<  qxxo .  25)  XENSUM)/(  1 . 0+DENSUM) 

Vd)  =  SqRT(d  .0-RKX0MEGAd)xx2)X(1.0-(OMEGAd)xx2/RK))) 

ADd)  =  1.0/(OMEGAd)XX2) 

DN=  1 . 0  +  (  SIGMAOXOMEGAd  )  )XX2 

B0(  I ):( (SIGMAOXVd  )  )xxz+(0nEGA(I)XH)xx2)/(DNXK2) 
Bld):(2.0xSIGnAOXV(l))/DN 
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18  CONTINUE 

H0=1.0 

IFdN.EO.OnHEN 
DO  20  I-I,IR 

HO^HOXBOLD/AOCI) 

20  CONTINUE 
HO:HOX(10.0«X(-0.05KASUBP)) 

ELSE 

DO  21  I-I,IR 

H0=H0XB0(I)/A0(1) 

21  CONTINUE 
H0=H0XSIGMA0 

ENDIF 

R=DBLE(SR) 

HO=DBLE(HO) 

A01-DBLE(A0(1).) 

B01=DBLE(BO(1)) 

B11-DBLE(B1(1)) 

IF( ( ORDER . E4 . 2 ) . OR . ( ORDER . E4 . 3) )THEN 
APROTd.Ol^HOXAOl 
APROTd,l)=O.0D+OO 
APR0Td,2)  =  H0 
APROTd,3)  =  O.0D+00 

if;order.eo.2)then 

APROT(2,0)=BO: 

APR0T(2,1)=B11 
APR0T(2,2)=1 .OD+OO 
ENDIF 

IF(0RDER.E«.3)THEN 

SO=DBLE(SIGMAO) 

APR0T(2,0)-S0XB01 

APR0T(2,1)=B01+(B11XS0) 

APR0TC2,2)=S0+B11 
APR0T(2,3)=1 .OD+OO 
ENDIF 
ENDIF 

IF((0RDER.EQ.4).0R.(0RDER.EQ.5))THEN 

A02=DBLE(AO<2)) 

B02=DBLE(BO(2)) 

B12=DBLE(B1(2)) 

APROrd,0)  =  HOXA01XA02 
APROTd, 11  =  0.00+00 
APR0Td,2)=H0X<A0I+A02) 

APROTd,3)  =  O.0D+OO 
APROTd,4)=HO 
APROTd,5)  =  0.0D+O0 
IF(0RDER.E9.4)THEN 
APROTC2,0)=B01XB02 
APROT(2,1)=IB11XB02)+(B01xb12) 
APROTt2,2)=B01+B02+(BllxBI2) 

APR0T(2,3)  =  BU  +  B12 
APR0T(2,4)=1 .OD+00 
ENDIF 

IFCORDER.EQ.SJTHEN 

SO=DBl£(SIGHAOJ 

APROT(2,0)=SOXB01XB02 

APR0TI2,  1)  =  (SOX((BUxB02)  +  <B01XB12)))  +  (B01XB02) 
APR0T(2,2)=(S0x(B01+B02+(B11xB12)))+(B11xB02)+(B01xB12> 
APROT(2,3)=(S0x(B11+B12))+B01+BO2+<B11xB12) 
APROT(2,4)=SO  +  BU  +  B12 
APROT<2,5)=1.0D+00 
ENDIF 
ENDIF 

IF((0R0ER,Eg.6) .OR. (ORDER. E0.7))THEN 
A02=DBLE(A0(2)) 

A03  =  DBLE(AO(3)) 

B02=DBLE(BO(2)) 

B03=DBLE(BO(3)) 

B12=DBLE(B1(2)) 

B13=DBLE(B1<3)) 

APROTd,0)=HOXA01XA02XAe3 
APROTd,I)  =  0.OD+OO 
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' 

APROT(1,2)<HOK((A02XA03)'»(A01M(A024^A05))) 

NEH07930 

APROT(l,3)<0.0D+00 

NEH07940 

V 

APR0T(1,A)>H0X(A01+A02«A03> 

NEW07950 

APROTd.SliQ.OD-fOO 

NEH07960 

APR0T(1,6>>:N0 

NEH07970 

IF( (ORDER. E9. 6). OR. (ORDER. Eq.7))THEN  ' 

NEH07980 

APROT(2,0)*B01KB02KB03 

NEH07990 

APROT(2,1)<(B11XB02»B03)+(B01«((B02«B13)+(B12»B03))) 

NEH06000 

APROT(2,Z)*(B02XB05)')’(B11K((B02kB13)'»(B12xB03>))+ 

NEH08010 

C  (B01«((B02+B03)4^(B12MB13))> 

NEH08020 

APR0T(2,3>:B01x(B12+B13)4(BllK((B02+B03>+(B12xB13>))'f 

NEH06030 

C  ((B02»B13)4(B12mB03)> 

NEH08040 

APROT(2,A)<B014(B11x(B124B13))4((B024B03)4(B12kB13)> 

NEH080S0 

APR0T(2,S)-B114B124B13 

NEH08060 

APROT(2,6)*1.0D+00 

NEW08070 

ENDIF 

NEWOBOBO 

1F(0RDER.E0.7)TH£N 

NEH08090 

SO-DBLE(SIGMAO) 

NEHOBIOO 

APROT(2,7)=1.0D+00 

NEW08110 

APROT(2,6)*(SOXAPROT(2,6))+APROT(2.5) 

NEH08120 

APROT(2,5)=(SOXAPROT(2,5))+APROT(2,«) 

NEW08130 

APR0T(2,A)=(S0XAPK0T(2,A))+APR0T(2.3) 

NEW08140 

APROT(2,3)=(SOXAPROT(2,3))+APROT(2,2) 

NEW081SO 

APROT(2,2)=(SOXAPROT(2,2))+APROT(2,1) 

NEW08160 

APRQT(2,1)=(SOXAPROT(2,1))+APROT(2,0) 

NEW08170 

APROT(2,0)=SOXAPROT(2,0) 

NEW0S180 

ENDIF 

NEH08190 

ENDIF 

NEW0S200 

ENDIF 

NEH08210 

C 

NORMALIZE  THE  ELLIPTIC  PROTOTYPES 

NEH0B220 

IF(TVPE.E0.5)THEN 

NEW08230 

DO  25  1=1,2 

NEW08240 

DO  26  J=0. ORDER 

NEWOB2SO 

APROTd,  J)=APROT(I,  J)/((DS(}RT(R))X»J) 

NEW08260 

26 

CONTINUE 

NEW08270 

25 

CONTINUE 

NEH0S280 

ENDIF 

NEW08290 

C 

APR0T(TVPE=6)  ALLOWS  THE  INPUT  OF  AMY  TYPE  OF  PROTOTYPE  FILTER 

NEW08300 

C 

BY  INPUTING  THE  COEFFICIENTS  OF  S 

NEW08  0 

1F(TYPE.EI}.6)THEN 

NEW08J. 

WRITE(6,X) 

NEHD8350 

WRITE(6,X)  'THIS  OPTION  ALLOWS  YOU  TO  USE  ANT  ANALOG  > 

NEW08340 

WRITE(6,X)  •  PROTOTYPE  OF  ORDER  1-10  ' 

NEW083S0 

WR1TE(6,X)  'XX  COEFFICIENTS  MUST  BE  IN  DOUBLE  PRECISION  ' 

NEW08360 

WRITE(6,X)  •  FORMAT  (0.1234567890400)  XX  • 

NEH08370 

WRITE<6,X) 

NEW083B0 

WRITE(6,X)  'INPUT  THE  NUMERATOR  COEFFICIENTS  FIRST* 

NEW08390 

DO  50  1=0, ORDER 

NEW0S4D0 

WRITE(6,X)  'WHAT  IS  THE  COEFFICIENT  OF  SXX'.I.'T’ 

NEW08410 

READ(5,X)  APR0T(1,I) 

NEW08420 

50 

CONTINUE  , 

NEW08430 

WRITE(6,X) 

NEW08440 

CALL  EXCMS  CCLRSCRN') 

NEWD8450 

WRITE(6,X)  'NOW  INPUT  THE  DENOMINATOR  COEFFICIENTS' 

NEW08460 

WRITE(6,X) 

NEW08470 

DO  51  J=0, ORDER 

NEW084S0 

WRITE(6,X)  'WHAT  IS  THE  COEFFICIENT  OF  Sxx',J,'T' 

NEW08490 

READ(5,X)  APR0T(2,J) 

NEW0S500 

51 

CONTINUE 

NEWOSSIO 

ENDIF 

NEW08520 

C 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXKXKXXXXKXXXXXKXXXXKXXKXKKXXX 

NEW08530 

C 

XXXXXXXKXXXXXXXXXXXXKXXXXXXNKXXRXXXKKXXXXXXXXXXXXXKXXXKXXKXXKXXX 

NEW08S40 

C 

DEFINE  ELEMENTS  OF  THE  BILINEAR  TRANSFORMATION  MATRIX 

NEN08550 

C 

XXXXXXXXKXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXKKXXKXXXKXX 

NEH08560 

C 

FOR  ALL  THE  BILINEAR  TRANSFORATION  MATRICES,  THE  FIRST  COLUMN 

NEH08S70 

C 

HAS  ALTERNATING  4l'S  AND  -I'S.  THE  LAST  COLUMN  HAS  ALL  41'S 

NEW08580 

■■ 

DO  60  1=1, ORDER 

NEII08S90 

DO  61  J=0,I 

NEH0860a 

BILIHR(I,J,0)=(-1.0)XXJ 

NEH08610 

BILINRd,  J,I)<1 .0 

NEH08620 

61 

CONTINUE 

NEH08630 

- 

60 

CONTINUE 

NEH08640 

M 
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C 

C 

c 


c 


77 

76 

C 


78 

75 

C 

C 


81 

80 

C 

C 

C 


102 

101 

100 

C 

C 


C 

106 

105 

111 

112 


C 

C 

C 

C 

C 

C 


FIRST  ORDER  MATRIX  HAS  ALREADY  BEEN  LOADED  CONTINUE  INTERATION 
UNTIL  YOU  HAVE  THE  APPROPRIATE  MATRIX 
DO  75  K-2, ORDER 
KM1=K-1 

LOAD  ELEMENTS  IN  INTERIOR  COLUMNS,  ALL  RONS  EXCEPT  LAST  ONE 
DO  76  L=0,KM1 
DO  77  M=1,KM1 
MM1=M-1 

BILINR(K,L,M)>BILINR(KM1,L,M)+BILINR(KM1,L,HM1) 

CONTINUE 

CONTINUE 

NON  LOAD  THE  LAST  RON 
DO  78  N-1,KM1 

BILINR(K,K,Nl-BILINR(K,a,N))t(lC-1.0D+00)M(K'6N)) 

CONTINUE 

CONTINUE 

XXKXKliXXXXXXXXKKXXXXXNXXXKXXKXKXNKXXIIXXKKKXXKXXXKKXKKKXKXKHKKKKK 

THE  FOLLOWING  CLEARS  OPROT  PRIOR  TO  LOADING 
DO  80  K°l,2 
DO  81  L-0,10 
DPROT(K,L)=0.0 
CONTINUE 
CONTINUE 

KKXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXKXXXKKXKXXXXKXX 

DIGITIZE  THE  ANALOG  PROTOTYPE  USING  THE  BILINEAR  TRANSFORMATION 

XXXXXXXXXXXXKXXXKXXXXXXXXXKXXXXXKXXXXXXXXxXXKKKXKKKKKKXKKNXXKXXX 

DO  100  N=l,2 
DO  101  J=0, ORDER 
SUMR=0.0 

DO  102  1=0, ORDER 

TEMP=APROT(N,I)XBILINR(ORDER,I,J) 

SUMR=SUMR+TEMP 
CONTINUE 
DPROTLN. J)=SUMR 
CONTINUE 
CONTINUE 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXKXXXXXXXKKKXKXXXXXXXXKKKXXXXXXKX 

OPTION  TO  SHOW  DIGITAL  PROTOTYPE  COEFFICIENTS 
CALL  EXCMS  CCLRSCRN') 

WRITE(6,X)  'THE  DIGITAL  PROTOTYPE  HAS  BEEN  COMPUTED.  DO  YOU' 
WRITE(6,X)  'WISH  TO  CHECK  THE  PROTOTYPE  COEFFICIENTS! ’ 

WRITE<6,x)  '  XX  1-YES/2-N0  XX* 

READ(5,x)  look 
IFCLOOK.EQ.DTHEN 

MAKE  DENOMINATOR  COEFFICIENT  OF  ZXXORDER  =1.0 
DO  105  K=l,2 
DO  106  L=0, ORDER 

DPROT(K,L J=DPROT(K,L)/DPROT<2,ORDER) 

CONTINUE 

CONTINUE 

CALL  EXCMS  CCLRSCRN') 

WRITE(6,x)  'THE  DIGITAL  PROTOTYPE  NUMERATOR  COEFFICIENTS  ARE.' 

DO  111  N=0RDER,0,-1 
WRITE(6,90Z)  DPROTtl,N),N 
CONTINUE 
WRITE(6,X) 

WRITE<6,x)  'THE  DIGITAL  PROTOTYPE  DENOMINATOR  COEFFICIENTS  ARE*' 
DO  112  M=ORDER,0,-1 
WRITE(6,902)  DPR0T(2,M),M 
CONTINUE 
WRITE(6,X) 

WRITE<6,X)  ’XX  1  TO  CONTINUE  XX* 

READCS.X)  LOOK 
ENDIF 

XXXXXXXXXXXXKKXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXKXXXKXX 

NOW  DO  THE  DIGITAL  TO  DIGITAL  TRANSFORMATIONS 

XXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXKXXXXKXKXKXXKKXXXXXXXXXXXXXXXKX 

IFKPASS.  EO.  1  ).OR.  <PASS.E5.2))THEN 

DEFINE  ELEMENTS  OF  THE  LOWPASS  AND  HIGHPASS  DIGITAL-DIGITAL 
TRANSFORMATION  MATRIX 

DEFINE  THE  ELEMENTS  FOR  FIRST  AND  LAST  COLUMNS 


KEW086Sa 

NEW08660 

NEH08670 

NEH08680 

NEH08690 

NEW08700 

NEH08710 

NEH08720 

NEW087S0 

NEH087A0 

NEM08750 

NEW08760 

NEH08770 

NEH08780 

NEW08790 

NEH08800 

NEW08810 

NEW08820 

NEH08S50 

NEW08S4D 

NEH08850 

NEW08860 

NEW08870 

NEM08880 

NEW08890 

NEW08900 

NEW08910 

NEW08920 

NEWD89S0 

NEW08940 

NEH08950 

NEW0S960 

NEW08970 

NEW08980 

NEW08990 

NEW090D0 

NEW09010 

NEW09020 

NEWD9DS0 

NEW09040 

NEW09050 

NEH09D60 

NEW09070 

NGMOVOBO 

NEW09D9D 

NEH09100 

NEW09110 

NEW09120 

NEW091S0 

MEW091<i0 

NEW09150 

NEW09160 

NEW09170 

NEW09180 

NEW09190 

NEW09200 

NEW09210 

NEW09220 

NEW09230 

NEW092<i0 

NEW09250 

NEWD9260 

NEH09270 

NEH09280 

NEH09290 

NEHO9S00 

NEH09S10 

NEM09S20 

NEH09SS0 

NEH09540 

NEH09S5C 

NEH09360 
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FILE.  DFCADO  FORTRAN  A 


DO  170  1-1,10 
DO  171  J^O.I 

LPHPTR(I,J,0)«(((-1.0)XALPHA)I«»J) 

IMJ=I-J 

LPHPTR(I,IMJ.I)«LPHPTRLI.J,0) 

171  CONTINUE 
170  CONTINUE 

C  FIRST  ORDER  MATRIX  HAS  IFFN  DEFINED  ABOVE 
C  SECOND  TO  TENTH  ORDER  BELOH 
DO  175  K:2,10 
KMlsK-1 

DO  170  L-0,KM1 
DO  177  M=1,KM1 
MM1=M-1 

LPHPTR<K,L,M)*LPHPTR{KM1,L.M)+((IPHPTR(KM1, I, MM17 )»{-!. OD+00) 
CXALPHA) 

177  CONTINUE 

176  CONTINUE 

DO  178  N=1,KM1 
KMN=K-N 

LPHPTR(K,K.N)=LPHPTR(K,0,KMN) 

178  CONTINUE 
175  CONTINUE 

C  XXXXXXKXXKKXKXXXXXXKKXXXXXKXXKKKXXKKXKKKKXXKKIIXKlIXXXXXKXXXXXKXIflf 

C  CALCULATE  THE  FINAL  FILTER 

C  FOR  HIGHPASS,  THE  PROTOTYPE  MUST  SUBSTITUTE  (-Z)  FOR  CZ) 
IF(PASS.E(J.2)THEN 
DO  180  K^l.Z 
DO  181  L=0, ORDER 

DPR0TCK,L)=DPR0T(K,L)X((-1 .OIXXL) 

181  CONTINUE 

180  CONTINUE 

ENDIF 

DO  185  K=l,2 
DO  186  L’O.ORDER 
SUMR=0.0D+00 
DO  187  M=0, ORDER 

TEMP=DPROT(K,M)XLPHPTR<ORDER,M,L) 

SUMR=SUMR+TEMP 
187  CONTINUE 

DFLTR(K,L)=SUMR 
186  CONTINUE 

185  CONTINUE 
ENDIF 

C  XXXXXXXXXXXXXKXXXXXXXXXXXXXXXXKXXXXXXXXXKXXXXXXKXXXXXXXXXXXXKXXX 

C  BANDPASS  AND  BANDSTOP  FORMULAS 

I F< ( PASS . E9 . 3 ) . OR . ( PASS . E5 . A ) ITHEN 
C  LOAD  THE  BANDPASS/BANDSTOP  TRANSFORMATION  MATRIX 

C  FIRST  LOAD  THE  FIRST  AND  LAST  COLUMNS 

DO  225  1=1,10 
DO  226  J=0,I 
BPBSTRCI, J,0)=CXXJ 
12=2X1 
IMJ=I-J 

BPBSTR(I,1Mj, I2)=BPBSTR(I, J,0) 

226  CONTINUE 
225  CONTINUE 

C  LOAD  THE  REST  OF  THE  FIRST  ORDER  MATRIX 
BPBSTRCI, 0,1)=(-1.0D+00)XB 
BPBSTRC1,1,1)=(-1 .0D+00)XB 

C  LOAD  THE  2ND-10TH  ORDER  TRANSFORMATION  MATRICES 
DO  230  K=2,10 
K2=2XK 
K2M1=K2-1 
KH1=K-1 
K2M2=K2-2 
KM12=2XKM1 
XM12M1=KM12-1 

C  LOAD  THE  SECOND  AND  SECOND  TO  LAST  COLUMNS 
DO  231  L=0,KM1 

BPBSTR(K,L,l)=BPBSTR<KMl,L,l)*(BPBSTRtKMl,L,0)xt-X.0D+00)XB) 

BPBSTR(K,0,K2M1)=(BPBSTR(KM1,0,KM12)x(-1.0D+00)XB)+ 


NEH09370 

NEM09380 

NEH09390 

NEH09A00 

NEH09410 

NEH09420 

NEH09430 

NEH09440 

NEH09450 

NEH09460 

NEH09470 

NEH09480 

NEH09490 

NEH09500 

NEH09510 

NEH09520 

NEH09550 

NEH09540 

NEH09550 

NEH09560 

NEH09570 

NEH095S0 

NEM09590 

NEH09600 

NEW09610 

NEH09620 

NEH09630 

NEN09640 

NEH09650 

NEH09660 

NEMD9670 

NEHD9680 

NEM09690 

NEH09700 

NEN09710 

NEN09720 

NEW09730 

NEH09740 

NEH09750 

NEMOOTBO 

NEW09770 

NEW09780 

NEH09790 

NEH09800 

NEH09810 

NEH09820 

NEH09830 

NEW09840 

NEHD9850 

NEH09860 

NEH09870 

NEW09880 

NEI’ID9890 

NEW09900 

NEH09910 

NEW09920 

NEM09930 

NEM09940 

NEW09950 

NEIt09960 

NEH09970 

NEH09980 

NEM09990 

NEHIOOOO 

NEHlOOlO 

NEH10D20 

NEH10030 

NEH10040 

NEH10050 

NEH10060 

NEH10070 

NEH10080 
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FILE 

OFCADO  FORTRAN  A 

C  (BPBSTR(KM1.0,KH12N1)KC) 

NEW10090 

KML=K-l 

NEWIOIOO 

BPBSTR(K,KML,K2Ml)=BPBSTR(K,L.l) 

NEWlOllO 

251 

CONTINUE 

NEH10120 

C 

NOW  LOAD  INTERIOR  ELEMENTS  EXCEPT  FOR  LAST  ROW 

NEWlOlSO 

DO  256  L°0,KM1 

NEW10I60 

\ 

no  237  M=2,K2M2 

NEH10150 

NEH10160 

MM2=M-2 

NEH10170 

BPBSTR(K,L,M)=BPBSTR(KM1,L.M)+(BPBSTR(KM1,L,MM1)«(-1.0D+00)K 

NEW10180 

C  B)+(BPBSTR<KM1,L.MM2)KC) 

NEW10190 

257 

CONTINUE 

NEM10200 

256 

CONTINUE 

NEW1D210 

C 

NOW  FILL  THE  LAST  ROW 

NEH10220 

DO  258  N-1,K2M1 

NEH10250 

K2MN=K2-N 

NEW10260 

BPBSTRC^.F  r''?MN)  =  BPBSTR(K,0,N) 

NEW10250 

258 

CONTINUE 

NEH10260 

250 

CONTINUE 

NEWI0Z70 

C 

lIXXXXKKKItKIiXXXXKXkXXKXKKKXKMKXKXXKKKXKKKNKICXKKKXXXKKXRKKXXKKXXXK 

NEW10280 

C 

CALCULATE  THE  FINAL  FIl TFR- 

NEW10290 

C 

FOR  THE  BANDPASS  CASE  REPLACE  (Z)  WITH  (-Z) 

NEM10500 

IF(PASS.E0.5)THEN 

NEH10510 

DO  260  M=l,2 

NEW10520 

DO  261  N=0, ORDER 

NEW10550 

DPR0T(H,N)=DPR0T(M,N)X(<-1 .0D+00)XXN) 

NEW1036L 

261 

CONTINUE 

NEW10550 

260 

CONTINUE 

NEW10360 

ENDIF 

NEW10570 

ORD£R2=2XORDER 

NEW10580 

DO  265  K=l,2 

NEW10590 

DO  266  L=D,0RD£R2 

NEW10600 

SUMR=O.OD+O0 

NEW10610 

DO  267  M=0, ORDER 

NEW10620 

TEMP=DPROT(K,M)xBPBSTR(ORDER,M,L) 

NEW10650 

SUMR=SUMR+TEMP 

NEW10660 

267 

CONTINUE 

NEW10650 

DFLTR(K,L)=SUMR 

NEW10660 

) 

266 

CONTINUE 

NEW10670 

265 

CONTINUE 

NEW106S0 

ENDIF 

NEW10690 

C 

xxxxxxxxxxxxxxxx 

NEW10500 

c 

OUTPUT  OPTIONS 

NEW10510 

c 

xxxxxxxxxxxxxxxx 

NEH10520 

IF((PASS.E9.3) .0R.(PASS.Eg.6))THEH 

NEW10550 

END=2X0RDER 

NEW10560 

ENDIF 

NEW10550 

IFKPASS.EQ.D.OR.  (PASS.EQ.2))THEN 

NEWI0S60 

END=ORDER 

NEH10570 

ENDIF 

NEW10580 

c 

MANE  DENOMINATOR  COEFFICIENT  OF  ZXXEND  E9UAL  TO  1 

NEW10590 

DO  250  M=l,2 

HEW10600 

DO  251  N=0,END 

NEW10610 

DFLTR(M,N)=DFLTR(M,N)/DFLTR(2,END) 

NEW10620 

251 

CONTINUE 

NEW10650 

250 

CONTINUE 

NEW1D660 

CALL  EXCMS  ( 'CLRSCRN' ) 

NEW10650 

WRITE(6,X) 

NEW10660 

WRITE(6.x)  'THE  FINAL  DIGITAL  FILTER  ISt' 

NEW10670 

WRITEt6,x) 

NEI-110680 

WRITE(6,X)  ■  NUMERATOR  COEFFICIENTS! • 

NEW1D690 

DO  500  I=END,0,-1 

NEW10700 

(•IRITE(6,902)  DFLTR<1,I),I 

NEW10710 

500 

CONTINUE 

NEW10720 

WRITE(6,X)  '  DENOMINATOR  COEFFICIENTS.’ 

NEW10750 

DO  501  J=END,0,-1 

NEW10760 

WRITE(6,902)  DFLTR(2,J),J 

NEH10750 

501 

CONTINUE 

NEH10760 

HRITE(6,X) 

NEW10770 

WRITE(6,X)  •  XX  ENTER  1  TO  CONTINUE  XX’ 

NEW107S0 

READ(5,X)  LOOK 

NEW10790 

• 

- 

C 

XXXXXXXXXXKXXXXXXXXXXKXKXXXXXXXXXXXXXKXKXKXXXXXXXXKXXXXKXKKXKXXX 

NEW10800 

II 
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FILEt  DFCADD  FORTRAN  A 


C 


310 


320 


OPTION  TO  SPOT  CHECK  FILTER  OUTPUT  MAGNITUDE  AND  FREOUENCY 
CALL  EXCMS  CCLRSCRN') 

WRITE(6,*)  'THE  TRANSFER  FUNCTION  FOR  THIS  FILTER  HAS  BEEN' 
HRITE(6.*)  'COMPUTED.  DO  YOU  HANT  TO  SPOT  CHECK  THE  MAGNITUDE' 
HRITE(6,«)  'AND  PHASE  OF  THE  FILTER  OUTPUT  FOR  A  GIVEN  ' 
WRITE(«,K)  'FREQUENCY?' 

WRITE(6,K)  '  1-YES/2-N0  »»  • 

READ(5,K)  LOOK 
IFCLOOK.EQ.DTHEN 
CALL  EXCMS  CCLRSCRN') 

IF(AFORDF.EQ.l)THEN 

HRITECO.X)  'HHAT  ANALOG  FREQUENCYCIN  KHZ)' 

WRITE(6,k)  '  DO  YOU  WANT  TO  CHECK? • 

WRITECO.K)  '  KK  DECIMAL  KX' 

READ(5,K)  F 
DF=(F/SSFREQ)K2.0«SPI 
ENDIF 

IF(AF0RDF.EQ.2)THEN 

WRITE(6,K)  'WHAT  DIGITAL  FREQUENCY  DO  YOU  WANT  TO  CHECK?' 
WRITE(«,«)  '  X*  DECIMAL  XX  • 

READ(5,X)  DF 
ENDIF 
X=COS(DF) 

Y=SIN(DF) 

Z=CMPLX(X,Y) 

XN=REAL(DFLTR(1,0)) 

YN=0.0 

XD=REAL(DFLTR(2,0)) 

YD=0 .0 

HOFZN=CMPLX(XN,YN) 

HOFZD=CMPLX(XD,YD) 

DO  320  N=1,END 
ZN  =  REAL(DFLTR(1,N))x(2)<»N) 

ZD=REAL(DFITR(2,H))x(Zxxn) 

H0FZN=H0FZN+2N 

HOFZD=HOFZD+ZD 

CONTINUE 

H0FZ=H0F2N/H0FZD 

FMAG  =  S()RT((REAL(HOFZ)»X2)  +  (AIMAG(HOFZ)XX2)) 

FMAGDB=20 . OXALOGIOC FMAG) 

FPHSE=ATAN(AIMAG<H0FZ)/REAL(H0FZ)) 

WRITE(6,*)  'THE  FILTER  OUTPUT  WILL  BE.' 

HRITE(6,91A)  FMAGDB 
WRITE(6,915)  FPHSE 
WRITE! 6, X) 

WRITE<6,*)  'DO  YOU  WANT  TO  CHECK  ANOTHER  FREQUENCY?' 

WRITE(6,X)  '  XX  1-YES/2-N0  XX  • 

READ(5,X)  look 
IFCLOOK.EQ.DGO  TO  310 
ENDIF 

KXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKKKXKXKXXXKKXXX 


OPTION  TO  LOAD  FREQUENCY  RESPONSE  INTO  A  DATA  FILE 
CALL  EXCMS  CCLRSCRN') 

'DO  YOU  WANT  THE  FILTERS  FREQUENCY  RESPONSE  ' 

'DATA  TO  BE  ENTERED  INTO  A  DATA  FILE  FOR  USE' 

'BY  A  PLOTTING  PROGRAM?  ' 

NOTE.  THE  FILE  NAME  WILL  BE  FILTER,  FILE  TYPE' 

WILL  BE  DATA,  100  POINTS  WILL  BE  CALCULATED' 

FOR  A  DIGITAL  FREQUENCY  FROM  0  TO  PI,  THE  FIRST' 
COLUMN  WILL  BE  THE  DIGITAL  FREQUENCY  (EXPRESSED 
AS  A  FRACTION  OF  THE  SAMPLING  FREQUENCY  0-1/2)  ' 
THE  SECOND  COLUMN  WILL  BE  THE  MAGNITUDE,  AND' 

THE  THIRD  COLUMN  WILL  BE  THE  PHASE! DEGREES) ' 

XX  1-MAKE  DATA  FILE/2-D0N"T  BOTHER  XX' 

LOOK 


WRITE!6,X) 

WRITE!6,x) 

HRITE(6,X) 

WRITE!6,X) 

WRITE!6,X) 

WRITE!6,X) 

WRITE!6,X) 

WRITE!6,X) 

WRITE!6,X) 

WRITE!6,X) 

WRITE!B,X) 

READ!5,X) 
IF!LOOK.EQ.l)THEN 
DO  350  1=0,100 
DF=REAL!I)X!5. 0/1000.0) 
DFR=REAL!I)X(SPI/100.0) 
X=COS!DFR) 

Y=SIN!DFR) 

Z=CMPLX!X,Y) 


NEW10810 
NEW10820 
NEW10830 
NEW10840 
NEW108S0 
NEW10860 
NEW10870 
NEH10880 
NEW10890 
NEW10900 
NEW10910 
NEH10920 
NEH10930 
NEW109A? 
NEH10950 
NEH10960 
NEW10970 
NEW10980 
NEW10990 
NEWllOOO 
NEWllOlO 
NEW11020 
NEW11030 
NEW11040 
NEWllOSO 
NEW11060 
NEW11070 
NEW11080 
NEW11090 
NEWlllOO 
NEWllllO 
NEW11120 
NEW11130 
NEW11140 
NEW11150 
NEW11160 
NEW11170 
HEW11180 
NEW11190 
NEW11200 
NEW11210 
NEW11220 
NEW11230 
NEWllZAO 
NEW11250 
NEW11260 
NEWn270 
NEW11280 
NEW11290 
NEW11300 
NEW11310 
NEW113Z0 
NEW11330 
NEWlllAO 
NEW11350 
NEN11360 
NEII1137  0 
NEW11380 
NEW11390 
'NEW11400 
NEW11410 
NEWllAEO 
NEWIIAIO 
.  NEWllAAO 
NEH11450 
NEWllASO 
NEWU<470 
NEW11480 
NEWU490 
NEW11500 
NEW11510 
NEH11520 
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FILE. 


351 


350 

901 

902 

903 
909 
905 
919 

915 

916 


DFCADD  FORTRAN  A 


XN=REAL(DFLTR(1.0)) 

YN=0.0 

XD=REAL(DFLTR(2,0)) 

YD=0.0 

HOFZN=CHPLX(XN,YM) 

HOFZD'CMPLXCXD.YD) 

DO  351  N=1.END 
ZN  =  REAHDFLTR(1,N))*(Z«»N) 

ZD=REAI(DFLTR(Z,N))K(Z«*N1 

HOFZN=HOFZN+ZN 

HOFZD=HOFZD+ZD 

CONTINUE 

HOFZ’HOFZN-'HOFZD 

FnAG2SQRT((REAL(HOFZ)KK2)+(AIMAO(HOFZ))»2>) 

IFCREALtHOFZ) .NE.O.OITHEN 
FPHSER=ATAN(AIMAG(HOFZVREAL(HOFZ)) 

ELSEIF( REAL (HOFZ) .EQ. 0.0 ITHEN 
IFIAIMAGf HOFZ) .GT. 0 . OlFPHSER'l .570796 
IFC AIMAG( HOFZ) .LT. 0.0) FPHSER=-1 .570796 
IF(AIMAG(HOFZ) .E8.0.0)FPHSER=0.0 
ENDIF 

FPHSE=FPHS£R»(iaO . 0/SPI ) 

NRITE(«,916)  DF.FMAG.FPHSE 
CONTINUE 
ENDIF 

■ixxirvKKXKKKXXKXXiCXXKXXKXKirxXXXXXXXXXXXXKKXKXKXKXXXXXXXilXXXXXNXXX 
FORMAT(5X,F10.5, •  SXX',11) 

FORMAT(5X,FI0.5, '  ZXX',11) 

F0RMAT(3X, 'ALPHAS  ',£9.9) 

F0RMATt3X, •<=  *,£9.9) 

F0RMAT(3X, 'THE  LOHPASS  DIGITAL  PROTOTYPE  FREQ.  IS  ’,£5.3) 
F0RMAT(3X, 'MAGNITUDE  =  ',F8.2,'  DB') 

F0RtlATt3X, 'PHASE  =  ',F6.2,'  RADIANS') 

FORMAT(3X,F9.3,10X,F5.3,10X,F6.2) 

STOP 

END 
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